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U ■ ABSTRACT 

We present results from a set of over 300 pseudospectral simulations of at- 
mospheric circulation on extrasolar giant planets with circular orbits. The sim- 
ulations are of high enough resolution (up to 341 total and sectoral modes) to 
resolve small-scale eddies and waves, required for reasonable physical accuracy. 
In this work, we focus on the global circulation pattern that emerges in a shallow, 
"equivalent-barotropic" , turbulent atmosphere on both tidally synchronized and 
unsynchronized planets. A full exploration of the large physical and numerical 
parameter-space is performed to identify robust features of the circulation. For 
some validation, the model is first applied to Solar System giant planets. For 
extrasolar giant planets with physical parameters similar to HD 209458 b — a pre- 
sumably synchronized extrasolar giant planet, representative in many dynamical 



-2- 



respects — the circulation is characterized by the following features: 1) a coherent 
polar vortex that revolves around the pole in each hemisphere; 2) a low number — 
typically two or three — of slowly-varying, broad zonal (east-west) jets that form 
when the maximum jet speed is comparable to, or somewhat stronger than, those 
observed on the planets in the Solar System; and, 3) motion-associated temper- 
ature field, whose detectability and variability depend on the strength of the net 
heating rate and the global root mean square wind speed in the atmosphere. In 
many ways, the global circulation is Earth-like, rather than Jupiter-like. How- 
ever, if extrasolar giant planets rotate faster and are not close-in (therefore not 
synchronized), their circulations become more Jupiter-like, for Jupiter-like rota- 
tion rates. 

Subject headings: planetary systems - planets and satellites: general - stars: 
atmospheres - turbulence 



1. Introduction 

A large number of extrasolar planets orbiting nearby sun-like stars is now knownQ So 
far, all of the planets are giant planets, which orbit ^ 5 AU from their host stars. These 
extrasolar giant planets possess very interesting and unexpected orbital properties. Some 
of them orbit extremely close to their host stars (semi-major axis a ^ 0.1 AU) and have 
very nearly circular orbits (eccentricity e ^ 0.05). These are the close-in extrasolar giant 



planets. If they possess tidal Q parameter values similar to that of . 
are expected to be in 1:1 spin-orbit resonance, or synchronous, states ( 


upiter (~10^), they 
Goldreich and Soter 


1966; 


Rasio et al. 


1996: 


Lubow. Tout. & Livio 


1997 




Oeilvie & Lin 


2004 


). The rest of the 



extrasolar giant planets orbit further away from their host stars, generally with substantial 
eccentricities. Among them, some have nearly circular orbits however, like the Solar System 
giant planets. In this paper, we focus our study on the atmospheric circulation of all giant 
planets with circular orbits — from the close-in extrasolar giant planets to the more distant 
"unsynchronized" extrasolar giant planets, as well as the Solar System giant planets. 



The existence of a dynamic, stably-stratified (radiative) outer layer is a common fea- 
ture on planets with atmospheres. Extrasolar giant planets are also expected to pos- 
sess such layers — independent of the proximity to their host stars. The dynamic layer is 



^See, e.g., \protect!ittp://exoplajiets.org/almcinacfrciine.html 

\protect||ittp : //www . obspm . f r/ encycl/encycl . html 



and 
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typically ^ 10 pressure scale heights thin, much smaller than the radius of the planet. 
It extends from just below the "weather layer" (at the bottom) up to the turbopause 
(at the top) — i.e., from near the top of the convection zone up to the tenuous region 
where molecular diffusion begins to dominate over advection. Following the first detec- 
tions of se condary echpses for three t ransiting extrasolar giant planets by the Spitzer Space 



Telescope (ICharbonneau et al.l l2005l : iDeming et al.l l2005l . l2006l ). the radiative layers of ex 
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Since vigorous motions and meteorology in these layers will strongly influence observable 
properties, a good understanding of atmospheric dynamics on extrasolar giant planets is 
needed. The understanding is especially crucial for the synchronized planets with their 
unique heating configuration. It is uncertain what temperature distribution should result 
from the combined ef fects of stationary dayside heating and advec tion on the atmospheres 
of these planets (e.g. IShowman and GuillotI l2002l : ICho et al.ll2003l ). This i s partly because 



they are dyriamica lly very different from the Solar System giant planets (jCho et al.l 12003 



Menou et al.l 120031 ). Close-in extrasolar giant planets are also independently important for 
theory since their expected synchronous states provide an idealized forcing situation, impor- 
tant for improving understanding of planetary (as well as possibly stellar) atmospheres in 
general. 

Recently, several global simulations of close-in extrasolar giant planet atmospheric cir- 
culation have been performed, focusing on di f ferent aspects of the overall complex problem 
(IShowman and GuillotI l2002l : ICho et al.l l2003l : ICooper and ShowmanI 1200 5l ) . Global calcula- 
tions are necessary to st udy the ful l effect of rotation on large-seal^ flows — the general 



circulation — on planets. ICho et al.l (120031 ) have carried out a high-resolution, one-layer 
turbulence simulation of the close-in extrasolar giant planet HD 209458 b using the pseu- 
dospectral method. They stress that the gross feature of the atmospheric flow on this 
"Hot Jupiter" is markedly different from that on the Jupiter: the horizontal circulation 
pattern is actually more like that of the Earth or Venus, ne ar and above the cloudtops. 
As elucidated in the companion paper by iMenou et al.l (120031 ), this is expected on dynam- 
ical grounds. The marked difference is due to the close-in planet's slower rotation rate 
(i7™ ^ 2.1x10"^ rad s~^) and hotter equilibrium temperature (T^ ^ r^ 1500 K), compared 
to tho se of Jupiter (1]J"p = 1.4x10-^ r ad s"^ and T.^^p = 130 K). IShowman and Guillot 



(120021 ) and ICooper and ShowmanI (120051 ) have also performed global simulations. Using the 
finite difference method in multiple-layer models, they focus on the laminar state and the 



^By "large-scale", it is meant horizontal scale L > Rp/IQ, where Rp is the radius of the planet. 
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vertical structure of the circulatio n. Both of these multiple-layer studies solve a more gen- 
eral set of equations than that of ICho et al.l (120031 ) . allowing vertical variations over many 
pressure scale heights to be represented. The multiple-layer studies also differ in their model 
setup: in most cases, an atmosphere at rest is heated, leaving out the important influences 
of pre-existing jet streams and evolving turbulent eddies. 

The structure and robustness of well-resolved horizontal flow forming under applied 
heating is the primary subject of this paper. Specifically, global circulation patterns possible 
in a stably-stratified, turbulent atmosphere on extrasolar giant planets with circular orbits 
are studied. This paper presents several much-needed developments in extrasolar planet 
characterization study: 1) a broad parameter-space exploration of a common type of ex- 
trasolar giant planet — i.e., extrasolar giant planets with circular orbits, which include the 
transiting close-in extrasolar giant planets; 2) a validation of the circulation model against 
a previously successful giant planet circulation model; and, 3) a detailed description of the 
emplo yed model so tha t equitable comparisons can be made with alternative flow rnodels . 
As in ICho et al.l (120031 ) . we here use the "equi yalent-barotro pic fo rmulat i on" ( Salbyl Il989l ) 
of the primitive equations of meteorology (cf., HoltonI ( 1992 ) and Salby ( 1996 )). hereafter 
referred to as the equivalent-barotropic equations. The equivalent-barotropic equations allow 
us to practically perform the large number of high-resolution simulations necessary for the 
exploration while bypassing the issue of the current lack of information on the vertical struc- 
ture (e.g., distribution of radiatively-active species). The equivalent-barotropic equations 
also allow a physically clear interpretation of heating and cooling associated with dynamics 
in a stratified layer. In this paper, we focus on the adiabatic dynamics — appropriate for 
when the dynamical timescales are shorter than the radiative timescales or when diabatic 
heating and cooling nearly cancel in the net. The focus sets the stage for future diabatic 
(non-adiabatic) calculations. 

The overall plan of the paper is as follows. In §2, we briefiy describe the current state 
of understanding of giant planet atmospheric circulation, laying stress on the aspects that 
are important for extrasolar giant planets. In §3, we present the primitive equations, the 
equations that govern the large-scale dynamics of atmospheres, and discuss the origins and 
properties of the reduced version of the primitive equations, the equivalent-barotropic equa- 
tions. Readers wishin g more details on the derivation of the equivalent-barotropic equations 
are referred to ISalbyl (119891 ). Readers familiar with the equations of atmospheric motion 
may skip the early part of this section. Later in the same section, we describe our repre- 
sentation of the differential thermal forcing, important for modeling close-in extrasolar giant 
planet atmospheres. In §4, we present the simulations of Solar System giant planets using 
the equivalent-barotropic equations. The simulations provide model validation as well as an 
initial assessment of circulation on unsynchronized extrasolar giant planets, which are essen- 
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tially modeled in the present work as Solar System giant planets with greater equilibrium 
temperature Tg and/or slower rotation rate fl. In §5, we summarize the detailed parameter- 
space exploration of the circulation on HD209458b, a representative close-in extrasolar giant 
planet in many dynamical respects. This brackets circulations possible on close-in extrasolar 
giant planets in the context of one-layer modeling. In §6, we conclude our presentation of 
the global adiabatic calculations. 



2. Atmospheric Circulation on Giant Planets 

The dynamically-active layer of a giant planet atmosphere is forced by agents both 
external and internal to the layer, which render the large-scale flow in the layer very complex. 
The forcing comes from i) the host star irradiation above, ii) the convection below, and in) 
the momentum and heat fluxes within. The impinging stellar radiation is mostly shortwave 
(i.e., ultraviolet and visible). Some of the incoming radiation is scattered by the atmospheric 
gases, which is composed mostly of inert H2 and He by volume. Some of it is reflected back 
by clouds, if they are present. And, some of it is absorbed by the atmosphere, particularly by 
active minor constituents such as H2O vapor and clouds. The atmospheric gases and clouds 
also emit and absorb longwave (i.e., infrared) radiation, leading to either further heat transfer 
between different parts of the atmosphere or heat loss to space. Convection provides direct 
stirring as well as latent heat from condensation. Roughly speaking, the static atmospheric 
structure is in radiative equilibrium above the weather layer (in the middle and upper part 
of the active layer) and in radiative-convective equilibrium below the weather layer (near 
the bottom of the active layer). In this way, radiative transfer processes are important 
in establishing the vertical temperature structure throughout the active layer. However, 
dynamical processes (such as transport and dissipation of eddied and waves) significantly 
modify the basic temperature structure that would be established by the radiative processes 
alone, in the absence of motion. 

On Solar System giant planets, despite the markedly different orbital, thermal, and 
chemical properties, the atmospheres all exhibit strong zonal (east-west) banding and asso- 
ciated zonal jets (longitudinally-averaged eastward winds), with long-lasting vortices em- 
bedded in them. The bands, jets, and vortices can all be understood as a dynamical 
equilibrium state of a st ably-stratified shallow layer of turbulent fluid on a rotating sphere 



( ICho and Polvanilll996al ). In a stratified-rotating fluid, the motion is predominantly layer- 



"Eddies" are the residual of the dynamical field after the subtraction of the mean. In geophysical fluid 
dynamics, "eddies" are often used interchangeably with "vortices" and even "waves". 
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wise and two-dimensional (2-D). As a consequence, the flow exhibits a simultaneous inverse 
cascad e of energy up to the large scales and forward cascade of enstrophjfl down to the small 
scales (IKraichnaru 1 1 96 71 : iLeithlllQGSl : lBatchelorlll969l ) . This is in marked c ontrast to the usual 
forward cascade of energy (only) in 3-D turbulence ( Kolmo gorovl 1 1 941a| jblfl) . In 2-D turbu- 
lence, th e inverse cascade r oughly corresponds to growth of eddies (vortices) by successive 
mergers (IMc Williams! 119841 ) . In an isotropic condition, the eddies ultimately grow to the 
size of the domain. On the sphere, however, a fundamental anisotropy due to the geometry 
is present, if the sphere is rotating: Coriolis acceleration provides a restoring force on the 
growing eddies in the meridional (north-south) direction while none in the zonal direction. 
In this anisotropic condition, eddies are free to grow to the largest available scale in t he zonal 
direction but restricted in the meridional direction — thus forming the bands and jets (IRhines 



19751 ). In this heuristic explanation of the bands, the characteristic meridional width of the 
bands is loosely defined by the Rhines scale, = 7r(2[///9)^/^, where U is the root mean 
square wind speed of the jets and P = 2Q cos (j)/Rp with cf) the latitude and Rp the radius of 
the planet. According to this mechanism, the number of jets and bands is ^irRp/Lp and is 
an increasing function of the planetary rotation rate Q. 



Cho and Polvanil (jl996bl ) have explicitly shown in numerical simulations that the banded 



structures on all Solar System giant planets can be reproduced very well through this 
mechanism". In their study, they use the shallow-water e quations, which formally applies 
to a thin layer of water with a free surface (jPedloskyl 119871 ) . The shallow- water simulations 
capture reasonably well the number, size, and amplitude of the zonal jets — the gross fea- 
tures of the atmosphere — on all four giant planets of t he Solar System. In additio n, the 
model captures the anticyclonic[f] dominance on Jupiter ( Mac Low and IngersoU 1986 ). It is 
important to note that their simulations use only the observed physical parameters values 
of fl, g, Rp, U, and He (respectively, planetary rotation rate, surface gravity, planetary ra- 
dius, global root mean square wind speed, and thickness based on the scale height at the 
cloudtop level). The qualitative match obtained between the simulations and the observa- 
tions gives confidence in the basic approach of using a barotropic, adiabatic, shallow-layer 
turbulence model to study giant planet atmospheric circulations. Similar results using the 
equivalent-barotropic equations is discussed in §4.1[ Alternative models of bands and jets on 
a giant planet also exist, however. One model is based on the surface expression of putative 
deep, convective, Taylor-Proudman columns aligned with the planetary rotation vector Q 



"'Enstrophy is the "vortical energy" of the fluid, defined by ^C^, where C = = Vxv is the vorticity 
(parallel to the vertical direction k) and v G is the horizontal velocity. 

^Cyclones (anticyclones) are vortices defined by the positive (negative) sign of C • where fl is the 
planetary rotation vector (and defines the north pole). 
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(cf., ISun et al.l (Il993l ). and references therein). Recently, iHeimpel et al.l (120051 ) have also 
obtained a good match to Jupiter's surface zonal wind structure using a 3-D global convec- 
tion model. Their calculation is performed in a thin spherical zone, which in principle is 
similar to the shallow layer turbulence approach (rather than the deep convection laminar 
approach) . 

Given the relative success of the shallow-water modej^, it is natural to try and extend 
the barotropic approach to extrasolar giant planets. However, a number of issues critical 
for the study of extrasolar giant planet atmospheric dynamics is currently either poorly 
understood or await observational constraints. This entails that current extrasolar planet 
atmosphere models, including radiative transfer models, necessarily contain some basic and 
crucial unknowns. The vertical distribution of active species has already been mentioned; 
Q is another. The la tter is crucial in characterizing the planet's visual, spectral, and hydro- 



dynamic propertie s (jCho and Polvanil Il996bl : IShowman and GuillotI l2002l : ICho et al.l 12003 



Menou et al.ll2003l ). Since orbital periods can be accurately measured, the rotation rate of 
close-in extrasolar giant planets can at least be inferred by assuming synchronization. How- 
ever, even in this case, there is still the question of how to initialize the model calculation. A 
background zonal flow, for example, will nontrivially affect the evolution of any developing 
flow by nonlinearly interacting with it. Without adequate initial and boundary conditions, 
it cannot be stressed enough that current model calculations should more properly be con- 
sidered as revealing important mechanisms and plausible flow states, rather than predicting 
reahty. 

With the above caveats in mind, we forge ahead by identifying giant planet properties 
that should be general. The flow in the radiative layer is generally not driven by the stellar 
heating directly but by the deviation from the equilibrium temperature. Moreover, the 
deviation responsible for the transport is caused by the atmospheric dynamics — specifically, 
by eddies and waves. Hence, it is more proper to think of the dynamics as controlling 
the temperature distribution, rather than the other way around. This is why it is crucial to 
represent the flow dynamics properly, both initially and throughout the flow evolution. Well- 
known examples of eddy/ wave-induced temperature and flow modifications exist. Rossby 
( "planetary") and gravity (buoyancy) waves cause major warmings over the cold pole, forcing 
the region to be above the radiative equilibrium temperature and undergo cooling. Waves can 
also induce quasi-periodic flow, and its associated temperature reversals, in the equatorial 
region. Cold regions at the subsolar point just above the troposphere and at the top of 
the radiative layer are eddy/ wave-driven as well, as the rising air in those regions cool 



^We must keep in mind, however, the gross simphfications inherent in the shallow- water model (see §3). 
The model does not — and in many cases cannot — address a number of important aspects of the circulation. 
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adiabatically. These are fairly well understood phenomena in the Ear th's atmosphere, where 
considerable observational data is available (cf.. [Andrews et al.l (Il987l )). Remarkably, similar 
heating and cooling mechanisms appear likely to operate on extrasolar giant planets, as 
we shall show later (§4 and §5). On rapidly-rotating and unsynchronized extrasolar giant 
planets, meridional transfer of heat is expected; and, on slowly-rotating or synchronized 
extrasolar giant planets, zonal transfer of heat is also expected, in addition to the meridional 
transfer. 

Under certain circumstances , Newtonian cooling a ppro ximati o n ma y be fruitfully used 
to represent the net heating (cf.. [Andrews et al.l (119871 ) and ISalbyl (Il996l )). The approxima- 
tion is a simple representation of radiative effects on large-scale dynamics — a relaxation of 
the temperature field to a specified "equilibrium" distribution, which in fact depends in a 
complicated way on the atmospheric constituents and their spatial distributions. It is not 
known a priori in general. The approximation is crude in that it ignores vertical motion 
{dp/dt), assumes maximum temperature perturbation T' small compared to the equilibrium 
temperature (i.e., T' ^ 1500 K for close-in extrasolar giant planets), approximates the verti- 
cal gradient of the transmission functi on to be ^l/ij/n , and requires the environment to be in 
local thermodynamic equilibrium (cf.. [Andrews et al.l (119871 )). It is widely used in theoretical 
studies of the Earth's atmospheric dynamics, but caution should be exercised when liberally 
extending its use to extrasolar giant planets. 

In this paper, we focus on the dynamical implications of a specified temperature dis- 
tribution containing a large, steady day-night component. In the context of the Newtonian 
cooling approximation, the relevance of this steady regime can be understood in terms of 
the relative magnitude of the advective timescale, Ta^L/U, and radiative timescale, ~ 
pCp{dT / dz*)^^ / (4crT^). Here, L, U, p, Cp, a, T, T, and z* are the characteristic length, char- 
acteristic speed, basic density, specific heat at constant pressure, Stefan-Boltzmann constant, 
equilibrium temperature, transmission function, and log-pressure height (= i/p In (pref/p), 
with Hp the pressure scale height and Pref a reference pressure), respectively. When r,. ^ Tq, 
advection of heat is much faster than radiative heat exchanges and the dynamics is essen- 
tially adiabatic. Based on preliminary diabatic calculations we have performed with explicit 
Newtonian cooling, the cases when > in our simulations are not qualitatively very 
different from the ^ case — unless Tr ^Ta- Hence, adiabatic simulations are relevant 
over a wide range of x^. In the opposite limit, when ^Tq, the flow dynamics equations 
cease to be physically connected to the fast radiative forcing. As we argue in §3.4, in this 
case it is reasonable to model this quasi-instantaneous radiative forcing limit by imposing a 
fixed temperature structure. A detailed study, which focuses on the relaxation of this model 
constraint, will be presented elsewhere. 



- 9 - 



3. Circulation Model 

3.1. The Primitive Equations 

The motion of an atmospheric fluid element is governed by the primitive equations. 
The equations are so called because they constitute the starting point for studying large- 
scale dynamics in the atmosphere. The presumption in their use is that the small scales, 
which are not resolvable either observationally or numerically, are parameterizable within 
the framework of large-scale dynamics. In atmospheric studies, p is commonly used as the 
vertical coordinate, rather than the geometric height, z = z{'x,p,t) with xg3?^; z becomes a 
dependent variable in this case. In the standard p-coordinates, the primitive equations are: 

Dv Mtan0, „ ^ „ ^ ^ /, X 

+ ^— — kxv = -Vp $ - /kxv + J^-P , (la) 



Dt R„ 



9$ 
dp 
duj 
dp 
DO 



Dt CpT 



-a , (lb) 



-Vp ■ V , (Ic) 



Qnet , (Id) 



where 



D d ^ d 

Dl = ^ + + 

In Eq. ([1]), v(x, t) = (m, f ) is the (eastward, northward) velocity in a frame rotating with O, 
and = gz is the geopotential, where g is the gravitational acceleration (assumed constant 
and includes the contribution from the centrifugal acceleration) and z is the height above a 
fiducial surface (e.g., the 1 bar pressure level); k is the unit vector normal to the planetary 
surface; / = 2f2sin(/) is the Coriolis parameter, the projection of the planetary vorticity 
vector 2Q, onto k; Vp is the horizontal gradient on a constant-p surface; uj = Dp/Dt is the 
vertical velocity; a is the specific volume, the reciprocal of density p; T and V represent 
momentum sources and sinks, respectively; 6 = T{pj.ei/pY is the potential temperatur^, 
where p^^i is the reference pressure and k = R/cp with R the specific gas constant; and, gnet 
is the net diabatic heating ratej^ The above equations are closed with the ideal gas law, 
p = pRT, as the equation of state. 



^The potential temperature is related to the entropy s hy ds — CpdlnO. 

®Note that various heating (>0) and cooling (<0) terms contribute to the diabatic heating rate q — e.g., 
stellar heating, longwave heating and cooling, latent heat, and sensible (conductive) heat. Newtonian cooling 
approximation is one simple parameterization of ^not- 
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To arrive at Eq. ([T]), one begins with the full Navier-Stokes equations — along with the 
energy equation , equation of state, and boundary conditions — in the spherical geometry 
(IBatchelorl 119771 ) . From the Navier-Stokes equations, two approximations are made at the 
outset . These are the "shallow atmosphere" and the "traditional" approximations (cf.. lSalby 
(Il996l )). The first assumes z/Rp -C 1 and is valid for the stable layers on all giant planets, 
including the extrasolar giant planets. The second is formally valid in the limit of strong 
stratification, when the Prandtl ratio [N/VlY ^ 1. In the ratio. 



N{z) 



9_ 
T 



(Tad-r) 



(2) 



is the Brunt-Vaisala (buoyancy) frequency, where Fad = g/cp and F = —dT{z)/dz are the 
adiabatic and the environmental lapse rates, respectively (Pedlosky 1987). The two ap- 
proximations allow horizontal (perpendicular to k) Coriolis terms to be dropped from the 
Navier-Stokes equations and vertical accelerations to be assumed small, as will be further 
described below. As an example, consider HD 209458 b: if we assume Jupiter-like com- 
position, temperature range of ~100 to ~1500 K across the globe, and nearly isothermal 
condition above the tropopausel^, we obtain (N/Q)^ ;^ 10^ with = g/{cpTf'^. The second 
approximation condition is also well met. Hence, as in many traditional atmospheric circu- 
lation studies, we may also reasonably take Eq. ([T]) as our starting point for extrasolar giant 
planets. 

As can be seen in Eq. (Pd), hydrostatic balance is assumed in the primitive equations. 
This condition restricts the vertical motion to be slow compared to the horizontal motion — 
or, equivalently, the vertical length scale of motions to be small compared to the horizontal 
length scale. It also allows local temperature to be associated with the atmospheric mass. 
That is, given two pressure surfaces (pi > P2) bounding a vertical region, the local thickness 
(or mass per area) of the r egion, = z^ — z^ , can be related to its average temperature via 
the hypsometric equation ([Andrews et al.l 119871 ): 



Az = -Hp{{T)) In ^ 



P2 



Pi 



(3) 



where Hp = R{T)/g is the local pressure scale-height and 



(T) 



dlnp 



P2 



-1 



Pi 



T(x, t) dlnp 



(4) 



P2 



^the level where the temperature ceases to decrease significantly with altitude (cf. JSeager and Sasselov 



(119981)) 
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is the region's average temperature. Note that Az, Hp, and (T) are all function of horizontal 
position and time. Eq. ([3]) is the basis of the pressure coordinate. Later, we will use an 
analog of Eq. to obtain horizontal temperature distributions on extrasolar giant planets. 



3.2. Equivalent-Barotropic Formulation 



For many applications, Eq. ([T]) is much too unwieldy and broad in scope. Moreover, in 
the absence of observational information on the atmospheric vertical structure to constrain 
the many parameters needed to i nitialize and so l ve Eg . ([H), a reduct i on of the equations 
is warranted — and necessary (cf., [Andrews et al.l (119871 ) and iPedloskyI (119871 )). A key ele- 
ment that must be represented faithfully when modeling the global circulation and thermal 
structure of planetary atmospheres is the dynamics of eddies and waves. Eddies and waves 
directly modify the flow and temperature structures on the large-scale by transporting mo- 
mentum, heat, and radiatively-active species. However, eddies and waves span in size from 
large scales down to scales well below those resolvable by most models solving the full primi- 
tive equations. A complete series of high-resolution 3-D calculations with all the key physics 
included and parameters explored is currently computationally prohibitive. Fortunately, 
some of this difficulty can be alleviated. In this work, we accomplish the reduction of the 
equations, resolving of the eddies and waves, and exploration of the parameter-space by 
verticall y integrating Eg. ([TD and using th e well-tested and highly-accurate pseudospectral 
method (IOrszag]|l970l : lEliassen et al.lll970l ) to solve the resulting, simpler set of equations. 
In doing so we are able to focus on some of the important elements critical for accurately 
capturing the global dynamics. 

There are other justifications for focusing on the column-integrated, horizontal dynamics 
as well. First, the region of extrasolar giant planet atmosphere that requires modeling for cur- 
rent observations is an extremely thin layer within the overall stratified layer, at most a few 
scale heights thick and located near the tropopause. It is in this sub-layer wher e most of the 
stellar irradiation is seen to be absorbed in radiative transfer calculations 



le.g. 



Seager et al. 



20051 : liro et al.ll2005l ) and from where the strongest (longwave) sensible features are likely to 
emerge obs ervationally. In addition, given the nearly isothermal vertic al structure in that 
region (e.g. ISeager and Sasselovlll998l : ISudarsky. Burrows fc Pintdl2000l ). we expect the fiow 
to be vertically well alignedj^"! In atmospheric dynamics studies, the use of reduced vertical- 
resolution models is common — especially for the tropopause region and above, where there 



^^Actually, good alignment may exist in the entirety of the stratified layer since the temperature structure 
is nearly isothermal throughout the layer. 
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is a strong b arotr o pic component (i.e. v ertical alignment; cf . . I Juckes and Mclntyrd ( 119871 ): 



Salby et al. ( 1990 ): Polvani et al. ( 1995 ): Ferreira and Schubertl (jl997 )). The typical length 
scale of the most unstable baroclinic mode (related to the non-alignment) is the Rossby 
defor mation radius, L-jj = NH/f ^ [gHy/"^ / f , where H Hp is the characteristic thick- 
ness ( jPedloskyl 119871 ). For baroclinic waves to exert a strong presence, L-ji should neither 
be -C Rp nor ^ Rp. As noted, the extreme smallness of L-ji has already been exploited to 
great advantage in mode ling Solar System giant planet circulations with a barotropic model 
(jCho and Polvanilll996bl ). For most unsynchronized extrasolar giant planets, we expect sim- 
ilar, or only slightly larger, values of L-ji {<^Rp). Even for close-in extrasolar giant planets, 
which receive intense heating, there is s ome justificatioii for using a barotro pic model since 
L-ji values for them are generally ^ Rp (ICho et al.ll2003l : iMenou et al.l 120031 ). In the future, 
it will be important to carefully delineate the regime of validity of our barotropic models 
with baroclinic (i.e., multi-layer) extensions. 



3.2.1. Equations 



Remarkably, ISalbyl (119891 ) has shown that both adiabatic (gnet = 0) and diabatic (g„ei 7^ 



0) forms of Eq. ([T]) directly reduce to a set of equations similar to the shallow-water equations. 



(Laplace 


1778; 


Gill 


1982) 



these dynamical equations reveal the clearest interpretation of the role of heating in isentropic 
coordinates, we have used that coordinate form of the equations in this work. In isentropic 
coordinates, constant ^-surfaces form the vertical coordinates, rather than p- or 2;-surfaces. 
The new, reduced set of equations — the equivalent-barotropic equations in 6'-coordinates — 
are obtained by assuming that the constant surfaces of thermodynamic variables {p, p, T, 9) 
share a common horizontal structure. A salient feature of the equivalent-barotropic model, 
as well as the shallow-water model, is that it is a global atmospheric dynamics model that can 
include the combined effects of differential rotation, stratification, and small-scale eddies and 
waves over long duration integrations. It also reproduces qualitative features of the global 
circulation on giant planets with a minimal set of adjustable parameters and assumptions. 
In this paper, we present results for the adiabatic case only, in which the effect of heating 
is parameterized by deflecting the modeled surface consistently with the temperature at the 
surface (see below and §3.3). A forthcoming companion paper will discuss the diabatic case. 

The equivalent-barotropic equations govern the dynamics of a semi-infinite gas layer, 
which is bounded below by a material surface. The bounding surface deforms according 
to the local temperature on the surface. In this work, the material surface is an isentropic 
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surface, 9o = constant. The governing equations read (ISalbylll989l ): 

Dv 



Dt 
DH 

'W 



where 



VIH + Hb] - /kxv + , (5a) 
kHV-v + Vn , (5b) 



D d ^ 

m = ^ + ^^^^ 



and 



The boundary condition, 



= f. (5e) 



D9o 

is trivially satisfied in the adiabatic case. In the diabatic case, this equation with forcing 
terms becomes a formal addition to the set of equivalent-barotropic equations. In Eq. ([5]), the 
dependent variables are the barotropic transformations of the original baroclinic variables in 
Eq. ([1]). For example, 

1 

v(x,)f:) = / w{yi,p,t) dp, (7) 

VO ^po(x,t) 

for the baroclinic velocity, v(x,p, t). Other variable have been similarly transformed. An 
equivalent-barotropic structure function, A = A{9{p)), is defined such that the inverse trans- 
formation for V is 



with the normalization. 



v(x,p,t) = A{e)v{x,t), (8a) 
rPo(x,t) 



1 

Po Jo 



A dp = 1. (8b) 



In Eqs. ([?]) and (IHaj) . barotropic and baroclinic forms are both represented by the same 
variable, but that should not elicit confusion. The P's in Eqs. (E]) represent dissipation, as 
before. The subscript "0" (as in po or 6q, for example) refers to the value of the variable at 
the bottom bounding surface. Accordingly, = -^(^o) is a measure of the local baroclinicity 
(i.e., vertical shear, or lack of vertical alignment), at the bounding surface. For an equivalent- 
barotropic structure which decays (grows) vertically, must be greater (less) than 1 to 
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preserve the normalization condition. In strongly-stratified atmospheric regions, in which 
the flow structures tend to be vertically aligned, should be close to unity. We have varied 
^0 from 0.8 to 1.2 and find no qualitative changes in the results^ In all the simulations 
presented in this paper, Aq = 1. In Eq. (jS^), zq = zo{-x,t) is the prescribed elevation of the 
bounding surface. Using the definition of 6', we have: 

Thus, the thickness, TC = TCi^Xjt), is essentially the temperature at the bounding surface. 
It is closely related to the local pressure scale height evaluated at the bounding surface, 
HpQ = RTq/q, and defines the potential temperature scale height, Hqq = H^qIk. It can be 
readily seen from Eqs. ([5]) that prescribed heating forces the fiow through the defiection of 
the bounding surface, which advects the temperature. The advected temperature in turn 
further drives the fiow, when gradients form. 

When «; = 1, the set of Eqs. (13) is fo rmally identical to the shallow- water equations, with 



bottom topography (cf.. lPedloskyl (119871 )). From the point of view of shallow- water dynamics, 
K (which is always < 1) primarily produces enhanced advection. As in the shallow- water 
equations, Eqs. ([5]) admit an important conservation law for the potential vorticity q: 

^ = '^1^ (IN 



where 

q (x, t) = 



(10b) 



and = VqiV^.V-j-i) is potential vorticity dissipation. In the absence of dissipation, the 
conservation law becomes: 

^ = 0- (11) 

Dt ^ ' 

Eqs. ffTOl) and ffTTl) are the equivalent-barotropic generalization of the shallow-water potential 

vorticity conservation. According to Eq. ( fTTl) . q is materially conserved in the absence of 

dissipation and can therefore serve as a proper tracer of the fiow. In this work we will 

make frequent use of q to visualize the fiow. In observed atmospheres, q is conserved very 

well, especially near the tropopause and above. Isentropic maps of q have been used very 

effectively to gain much understanding in atmospheric dynamics studies of the Earth (cf.. 



11 



In general, the results are sensitive to large differences in ^o- 
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Hoskins et al.l (Il985l )). Indeed, one of the significant advantages of our approach is our 



abihty to conserve this crucial quantity to a high degree. This is achieved through the high 
resolution of our calculations and the accuracy of the employed numerical method, which we 
now briefiy describe. 



3.2.2. Numerical Solutions 



Eqs. ([5]) are solved numerically in full spherical geometry — with x= (A,</)), where A is 
the longitude and is the latitude. The equations are solved in the vorticity- divergence 
form — i.e., the curl and divergence of the Eqs. (jS^), along with Eq. dSb)- In this form, the 
field variables, (C, 5, Ti) = (k- Vxv, V;v, HiTo)), are projected o nto the Legendre basis via 
the Gauss-Legendre transform (IOrszag||l97Cll : lEliassen et al.lll97Cll ): 



e(A,0,t) 



N=N(M) n 



im\ 



(12) 



n=l 



where ^ represents an arbitrary field variable, is the associated Legendre polynomial for 
total wavenumber n and sectoral wavenumber m, and N and M are the total and sectoral 
truncations. When = M, the truncation mask is a triangle in spectral space, and the 
resolution of a calculation is referred to as "T106", for example, if = 106. In physical 
space, this truncation correspo nds to a mini mum of (longitude, latitude) = (3A^ + 2 = 320, 
(3A^ + 2)/2 = 160) grid points flOrszadll970l ). In most spectral calculation is com- 

putationally much mo re accurate (though less efficient), per degree of freedom, than a finite 
difference calculation (jCanuto et al.lll988l : iHaltiner and Williamslll980l ). The high resolution 
and accuracy achieved in the spectral algorithm is crucial for allowing small-scale structures 
and turbulence to evolve without suffering much dissipation. A small hyperviscosity, of the 
form (— l)^"'"^z/V^^ for x = {1)2,3,4,8}, is included in Eq. ([5]) to stabilize the numerical 
integration while extending the turbulent inertia l range . Details on the use and effects of 
hyperviscosity can be found in ICho and Polvanil (jl996al ). As for the time-integration, the 
second-order accurate leapfrog scheme is used. We have thoroughly explored the numerical 
parameter-space associated with the present work and identified the region in which our 
results do not depend qualitatively on the precise choice of numerical parameter values. 

Throughout this paper, simulations described are initialized with a random turbulent 
fiow. This initialization crudely represents the effect of stirring arising from a variety of 
barotropic and baroclinic processes; its use is nearly universal in turbulence simulation stud- 
ies. The stirring is applied only at the start of the simulation and does not represent a 
continuous source of small scale forcing: it is used to generate a reasonably constrained 
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background flow in our atmospheric simulations. The motivation behind this key ingredient 
of our modehng strategy is further discussed in §3.4. The random vorticity field is character- 
ized by a narrow band in spectral space with a specified mean amplitude that corresponds 
to a global root mean square velocity U. The divergence field is set to be uniformly zero ini- 
tially. Given these fields, the remaining variable, Ti, is obtained by requiring that the initial 
flow be nonlinearly balanced — i.e., satisfy the divergence of Eq. ([5)d) in the limit d/dt 0. 
This procedure allows the turbulent calculation to proceed over long times without "blowing 



up" . A detailed discussion of the initialization procedure is also given in ICho and Polvani 



(jl996al ). We have verified that the results are insensitive to the location of the initial distri- 
bution's peak in spectral space, as long as the location is much smaller than the scale defined 
by Ljs (Rhines scale). The mean amplitude of the initial perturbations is a free parameter 
of the model, which we express here in terms of tJ. 

The equivalent-barotropic model admits the following parameters: ^o, ^, R, ^, 9, Rp, U, 
and To- Here, Tq is the global mean Tq. Note that three parameters appear in the equivalent- 
barotropic model, in addition to the ones that are present in the shallow-water model. 
They arise from the more explicit representations of the physics in the equivalent-barotropic 
model, compared to the shallow-water model. The added cost of the extra variables is 
minimal, however, since nearly all of the needed parameters can be reasonably constrained 
by observations or plausible scaling analysis. 

In the present work, if Jupiter is used as a paradigm giant planet to constrain R and 
Q (for the case of unsynchronized extrasolar giant planets), f7 becomes the most uncertain 
of the parameters. For Jupiter, R is 3779 J Kg~^ K~^, using k = 2/7; we have found that 
a K value of 2/5, appropriate for monatomic gases, also does not qualitatively change the 
results. In addition, f/ ^ 70 m s^^ for Jupiter; note that U in the upper tropospheres of all 
the Solar System planets is known and is between ~30 m s~^ and ~400 m s~^. However, the 
appropriate value for current extrasolar giant planets is not known and not easily constrained, 
leading us to adopt a somewhat-uncomfortable wide range: ~10^ m s~^ to ~10^ m s~^. While 
not stultifying, the uncertainty is not insignificant. As we shall show, values in the lower 
range can lead to identifiably different circulation pattern than those in the upper range. As 
yet, we are unable to definitively conclude on the exact circulation pattern without further 
observational constraints. Currently, there is no way to obtain U from first principles. 



3.3. Thermal Forcing 



Two points are important to note when addressing thermal forcing on extrasolar giant 
planets: i) in general, it is not the direct stellar irradiation but the net flux (i.e. that which 



-ir- 



is non-linearly modified by the circulation) that forces the temperature tendency (Eq. ([Hi)) 
and balances the flow; and, ii) because some of the stellar irradiation will be rejected by 
the atmospheric gases and clouds, the radiative equilibrium temperature distribution is not 
known with certainty for extrasolar giant planets at the present time. Further observations 
are required. Note that ii) necessarily implies that the strength and distribution of forcing 
in i) is unknown for extrasolar giant planets and must be treated as a parameter. 

For the Solar System giant planets, for which good flux observations exist, it is rea- 
sonable to model their atmospheres as layers of uniform thickness (temperature). This is 
because the emergent flux (which is as much as, or larger than, the insolation) is nearly 
uniform on them and ^ r^. Indeed, the fiux-associated equator-pole temperature gra- 
dient on Jupiter is only ~3 K at all longitudes. This suggests that thermodynamics ma y 



be effectively decoupled from the dynamics, as has been done in ICho and Polvanil (Il996bl ). 
Such a prescription is also probably appropriate for many unsynchronized extrasolar giant 
planets, which may behave plausibly like a "warm Jupiter" . However, the uniform thickness 
assumption is not likely to be as valid for the extremely close-in giant planets, given the much 
smaller than for Jupiter and the expected high radiative-equilibrium temperature gradi- 
ents. On these planets, the true temperature at the substellar point should be closer to the 
equilibrium temperature than in the polar — and if synchronized, antistellar — regions. Hence, 
an estimate of the radiative-equilibrium temperature distribution is required to model the 
flow on close-in giant planets. 

For synchronized planets, in the absence of circulation, the radiative equilibrium tem- 
perature on the day side is modeled as 

Teq(A,0) = T.,(^^y\l-A,^f'icosXcos<pY/\ (13) 

where T^, and i?* are the parent star's effective temperature and radius, respectively. Here, 
D is the distance of the planet from the parent star; and, Ax^ is the local value of the 
albedo, a function of the location on the day side. This definition differs from the standard, 
"averaged" one, often found in the literature: 



1/2 



/(l-^) (14) 



where A is the geometric albedo (an averaged value over the entire day side) and / is a 
parameter that describes whether the stellar flux is re-emitted by the entire planet (= 1) or 
only the day side (=2). Note that because the albedo is a function of temperature and space 
(via cloud formation, for example), the radiative equilibrium temperature field itself is not 
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really well defined — even in the absence of motion. However, these two definitions allows us 
to very roughly bracket the range of forcing amplitude that needs to be explored. 

Prior to synchronization, the outgoing flux was more likely "homogeneously distributed" 
over the planet's surface. There is no orography or topography to make the diabatic heating 
zonally asymmetric on a gaseous giant planet. The temperature resulting from a planet which 
emits homogeneously is given by Teq in Eq. (fT4l) with f = l. After reaching synchronization, 
however, the incoming stellar flux distribution is not zonally-symmetric and the equilibrium 
temperature associated with this influx is not uniform. The maximum radiative equilibrium 
temperature is then given by Teq (A = 0, = 0) in Eq. f|T3l) at the substellar point with the 
temperature profile decreasing as (cos A cos (pY^^ away from the point. Statically, the thermal 
forcing resulting from synchronization can therefore be modeled in first approximation as a 
permanent temperature perturbation superimposed on a layer of nearly uniform temperature, 
Teq. If the differences in local vs. global albedo values are ignored (i.e., A = assumed) 
for the sake of argument, Eqs. (fT3|) and f|T^ indicate that the amplitude of thermal forcing 
at the substellar point could be as large as ~ 40% in excess of the uniform temperature 
value, Teq (obtained for / = 1). Given the above-mentioned uncertainties in the definition 
of radiative equilibrium itself, we will consider amplitudes of thermal forcing reaching up to 
20-40% in our circulation models. 

On the night side, absent any circulation, atmospheric temperature should be approxi- 
mately the effective temperature determined by the internal luminosity. The internal lumi- 
nosity is powered by the slow gravitational contraction and possible helium rainout in the 
plane© 



Teff ~ lOOK (tW^ ' (15) 

where Tint,jup ~ 8.7 x 10^^° J s^^ is the estimated internal luminosity of Jupiter. With the 
observed temperatur e at the substellar po int of HD 209458 b possibly as high as ~1700 K (cf.. 



IP 

Ui 



Deming et al.l (120051 ) : ISeager et al.l (120051 ) ). the temperature difference between the substellar 
and antistellar points on this extrasolar giant planet could be as high as ~1600 K in the 
absence of motion. Obviously such an enormous temperature difference would be quickly, and 
violently, neutralized by atmospheric motion as long as < r,,. It is crucial to note that, if 
f2 7^ 0, the readjusted flow state will not be that which merely balances the above temperature 
gradient but also one which is very stro ngly modi fied by the planetary rotation. This is a 



fundamental property of rotating fluids (jGilllll982l ). Note also that is locally a nonlinear 



function of T and of the motion. Venus, which has a remarkably homogeneous temperature 



1^ Other mechanisms may be at play as weh (cf., iGuillot et al.1 (|2006f )). 
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distribution despite its extremely slow longitudinal variation of insolation, may serve as 
an example of such a post- adjusted distribution. Incidentally, during such an readjustment 
phase, the flow state may not be unlike the turbulent initial condition used in our simulations, 
since small-scale gravity waves and turbulence will be generated during the process. 

The salient point here is that a reasonable response to irradiation is to make the atmo- 
spheric layer thicker (due to the warmer radiative equilibrium temperature) on the day side 
than on the night side, as dictated by Eq. Accordingly, in this work, the flow is forced 
with heating by appropriately "puffing up" the atmospheric layer on the heated side while 
maintaining a global equilibrium temperature Tq (= Tgq)- This is accomplished by intro- 
ducing a permanent deflection of the bounding material surface through a specifled zo{x.) in 
Eq. ([Sfe). In most runs, the deflection is "grown" after the basic flow pattern is established 
with a characteristic e-folding time (usually ~5r, where r = 27r/i7 is the rotation period) 
in order to minimize introducing a large unbalanced component (small-scale fast modes) in 
the flow. A sudden "turning-on" of the forcing introduces a large amount of fast modes, 
which are in general insigniflcant for met eorology and are flltered out in numerical weather 



predictions ( iHaltiner and Williams! 119801 ). No heat (thickness) is added overall, but merely 
redistributed by the flow. Generally, this is most appropriate when 3> or when gnet 
(Eq. (Pl)). Here, we are mainly concerned with the horizontally advected component of the 
temperature, neglecting exchanges with regions above and below the modeled layer. How- 
ever, as we argue in §3.4, because of the way in which the bounding surface is deflected, our 
steady forcing procedure may also capture important features of the opposite regime, in the 
fast radiative limit 0. 

The fluid layer is bulged on the day side according to: 

n{X,(f)) = n + 7^' cos cos A (16a) 

and 

Zo(A,0) = Zq — 2;qCos0cosA, (16b) 

where Ti.', Zq > are the constant, maximum perturbation amplitudes of forcing at the 
substellar point above the reference levels, H and zq. H' is set at the start of the simulation 
to ensure that the flow is down-gradient of the temperature distribution and then to 
thereafter. Without loss of generality, we set zq = 0, since it can be subsumed into H, but 
Zq is held flxed throughout the simulation. That is, a permanent deflection of the lower 
bounding surface is maintained while a mass flow initially away from the substellar point is 
ensured, consistent with heating. The same perturbation amplitude, but with opposite sign, 
is applied to the layer on the night side so that the forcing integrated over the entire planet 
surface 5* is zero and 

/ 7^(A,0) dcpdX = n. (17) 
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We have also performed simulations in which the day side merges smoothly into a flat night 
side and found that the results do not change qualitatively. According to the argument 
outlined above, if we define 



we are interested in rj < 0.4. Since the actual value of rj may be different for different extra- 
solar giant planets and different altitudes, we vary its value from to 0.4 in our simulations. 

Note that the deflection approach cannot be used in the much lower levels of extrasolar 
giant planets since the vertical entropy gradient becomes very shallow in the equatorial region 
and a 6'-surface cannot be practically defined there0 The bulging associated with heating 
can in fact be represented either indirectly by lowering the height of the lower bounding 
surface (thus producing a "surface deficit") or directly by adding mass (thus increasing H). 
Similarly, the reduced thickness of the fluid layer on the night side (due to cooling) can be 
represented by directly removing mass, or by increasing the height of the lower bounding 
surface. This is just a restatement of the hypsometric relation, Eq. ([3]). The direct method 
is "diabatic" and the indirect method is "adiabatic". We have investigated both methods 
but present only the adiabatic case in this paper. Extensive study of the diabatic case will 
be presented elsewhere. Our simulations also indicate that the key results on the developed 
circulation do not depend qualitatively on the precise functional form of the differential 
thermal forcing — e.g., Gaussian or "cos^''''" dependence, as opposed to the simple "cosine" 
dependence of Eq. ( !T6|l . 

Because the bulge associated with thermal heating can introduce a large unbalanced 
component in the flow, as described above, it can quickly lead to a numerical instability. 
Recall that the pre-existing part of the flow is nonlinearly balanced. Unless the flow is "sta- 
bilized" with draconian dissipative measures, which are unphysical, a gentle "ramp up" to 
the equilibrium profile is necessary. Even in this case, the destabilizing effect of imbalance 
when the amphtude is large cannot be completely avoided. We have tested both instanta- 
neous and gradual forcing schemes. In the instantaneous version, the fluid is fully bulged 
right at the beginning of the simulation. In the gradual version, forcing is slowly increased 
with time during the numerical simulation, up to the chosen t] value. Again, we have found 
that simulations with instantaneous or gradual forcings gave results that are qualitatively 
similar. Overall, we find that the developed flow is surprisingly robust under a wide variety 
of conditions. This is one of our principle conclusions from this study. 



^•^Computationally, the large deflection required causes holes to develop in the model layer, leading to 
blow-up. 



n 
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3.4. Key Assumptions and Potential Caveats 

At this point, it is useful to review the main assumptions in our modeling strategy, in 
order to remind the reader of some of the potential limitations of the results presented in the 
ensuing sections. We model atmospheric circulation on close-in extrasolar giant planets with 
a single, adiabatic, equivalent-barotropic layer, which is steadily forced on the large scale as 
described in §3.3. Except in specific prescribed amount of small-scale turbulence is 

applied one time at the start of the model simulation. This is to generate a background flow 
which is consistent with one under no thermal forcing and which interacts with one forming 
due to our steady thermal forcing. Our aim here is to study, in a highly idealized way, the 
interaction between the background and forced flows. 

Focusing on equivalent-barotropic flows allows high enough spatial resolution to simulta- 
neously address the important flow interaction and explore the parameter-space well enough 
to ensure robustness of the results. While we have provided some justifications for the 
relevance of equivalent-barotropic character of the modeled region of the atmosphere, baro- 
clinicity due to vertical variations is significant and must be included for additional physical 
realism. Hence, it is important to extend the present one-layer calculations to multi-layer 
calculations. Ideally, high resolution should also be employed in the vertical direction for the 
same reason given for the horizontal direction. Extra care should be given to the specification 
and effects of vertical boundary conditions as well. 

The assumption of adiabaticity is clearly a simplification and one that is not applicable 
to all regions of the atmosphere. However, it is valid in many regions, where pressure is high, 
temperature is low, or heating and cooling cancel in the net. In addition, it allows the effects 
of diabatic forcing to be put in lucid context. In particular, the present study is a first step 
towards assessing the crucial interplay between "adiabatic wind" (as parameterized by U) 
and "diabatic wind" (as parameterized by t]). This step is significant because the relationship 
between wind structure and thermal forcing on planetary atmospheres is complicated and 
not completely understood: currently, there is no theory that explains the magnitude of the 
observed wind speeds at cloudtop heights. Hence, the wind speed cannot be naively related 
to the amount of incident stellar flux. Table [H showing a non-trivial scaling between U 
and semi-major axis a for Solar System giant planets, exemplifies this point. In this respect, 
atmospheric simulations which attempt to predict wind speeds purely on the basis of applied 
thermal forcing are physically limited, even though such model setups may be conceptually 
appealing. 

One regime of circulation for which our adiabatic simulations may turn out to be useful 
is the short radiative cooling time limit, — ^ 0. This may appear counterintuitive at 
first since the usual assumption of adiabatic condition corresponds to the opposite limit. 
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Tr —>■ oo. Notice, however, that in the hmit ^ 0, a parcel of air relaxing according to 
Newtonian "cooling" adjusts instantaneously to the state of radiative equilibrium (unless the 
"cooling" is balanced by other effects, such as conduction). If <^ r^, the sound crossing 
timesscale {=Rp/cs), such a fast forcing of the thermodynamic equation (Eq. [T]i]) cannot 
be balanced by global-scale motions. In this limit, one might sensibly filter the fast process 
out of the dynamics, similar to the way in which sound waves are filtered out by imposing 
hydrostatic balance or gravity waves are filtered out by imposing quasi-geostrophy. Note 
that when Tr^Ta, the adiabatic approach is not justified and explicit diabatic calculations 
are necessary. 

Another major assumption in our modeling approach is that extrasolar planet atmo- 
spheres are turbulent, in the sense that the flow is dynamically active on many scales, and 
that turbulent eddies and waves play a significant role in establishing the observed flow. 
Turbulence, as signified by the presence of eddies and waves, in the atmosphere is produced 
by a variety of mechanisms; and, it is present on all Solar System planets — even if visible 
tracers like clouds are not there to spotlight it. In our simulation set-up, one may question 
the validity of using a small-scale turbulent forcing for close-in extrasolar giant planets: given 
the expected dominance of large-scale external irradiation, one might argue for the absence 
of small-scale stirring — particularly from obvious sources such as convection and baroclinic 
instability. Here, we emphasize that the small-scale forcing in our simulations is not used to 
literally model those processes but used to facilitate a dynamically constrained atmosphere 
that is not at rest. After the initial one-time stirring, the flow in our close-in extrasolar 
planet simulations is forced only with a steady, large-scale forcing. 

In summary, we recognize the difficulty of assessing the applicability of our assumptions 
on extrasolar giant planetary atmospheres in the absence of almost any direct observational 
constraints. Since key ingredients necessary to reliably model the atmospheric dynamics of 
close-in extrasolar giant planets are presently not determined, we have adopted an approach 
that emphasizes the following: 1) previous successful modeling strategy, 2) physical processes 
that are reasonably well understood on Solar System planets, and 3) a baseline that allows 
effects of additional physics, such as diabatic heating, to be clearly delineated. It is standard 
to employ such an approach initially in atmospheric studies. 
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4. Unsynchronized Giant Planets 

4.1. Solar System Giant Planets 

Given the limited amount of direct information on atmospheres for extrasolar giant 
planets, validation of theoretical mo dels is critical. The curre nt work is closely-related to the 
shallow-water equations studies by ICho and Polvani Jl996al lb[). Their findings are directly 



relevant and lends some credence to our work. Hence, we briefly summarize the salient 
features before presenting results with the equivalent-barotropic equations. 

For shallow-water equations in spherical geometry, two dimensionless numbers define the 
dynamics for a fixed radius: the Rossby number, Ro = U/{fL), and the Froude number, 
Fr = U / {gHy^"^. These numbers arise from the ratio of terms in the momentum equation of 
the shallow-water equations. The numbers also appear in the equivalent-barotropic equations 
and control the dynamics in the same way. Ro is a measure of relative vorticity compared 
to the planetary rotation (or, the rotation time compared to the advection time), and Fr 
is the shallow-water analog of the Mach number in compressible fiows0 A convenient and 
important third dimensionless number can be constructed from the above two: the Burger 
number, = R^/F^ = L-ji' / L"^ . In rotating-stratified fluids, L-ji acts as an e-folding length 
limiting vortex and jet interactions. In planetary applications, the (3 (Rhines scale) effect 
is necessary, but not sufficient, to produce long-lasting banded structures; L-ji must also be 
finite. That is, the fluid must have a free surface, or be 2D-compressible . Quantitatively, 



Ly, mu st be ^ -Rp/3 for a banded structure to be stable. In their studies, ICho and Polvani 



Jl996al lb[) use the simplest form of shallow-water equations, without external forcing or 



topography, in order to clearly delineate the intrinsic flow evolution from that due to forcing. 
Only a delta-function force in time is used to stir the fluid at the beginning of the simulation. 
In the cases in which a simple, continuous Markovian random forcing is used, the physical- 
space evolution is strongly influenced by the forcing and vortices do not merge robustly 



fICho et al.lll997f ). 



Because of the close relationship, shallow-water equations results can be used to vali- 
date the current model. We do this by performing simulations of Solar System giant planet 
atmospheres using the adiabatic equivalent-barotropic equations. This not only gives con- 
fidence in the present code, it also sets the context for more complex diabatic calculations. 
In addition, direct comparisons with shallow-water equations simulations allow the effects 
of additional physical parameters in equivalent-barotropic equations simulations to be care- 



^"^In the shallow-water equations, surface gravity waves are admitted, not sound waves, since the fluid is 
assumed homogeneous and incompressible in 3-D. 
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fuUy assessed. Since the primitive equations require even greater number of parameters, 
rendering model results sometimes difficult to interpret clearly, equivalent-barotropic equa- 
tions simulations provide a good benchmark for full primitive equations models. The use 
of shallow-water equations model to elucidate or verify primitive equations model results is 
a common procedure in atmospheric dynamics studies. Moreover, present calculations of 
Solar System giant planets are useful in themselves since they have never been simulated 
with the equivalent-barotropic equations. The calculations also have value in that Solar 
System giant planets still serve as paradigmatic giant planets in many ways; there must 
surely be extrasolar giant planets which are Solar System giant planet-like in their atmo- 
spheric properties. Finally, we note that the equivalent-barotropic equations are preferable 
to the shallow-water equations because they derive directly from the primitive equations; 
certain parameters (e.g., L-ji) acquire a more physical, less arbitrary, interpretation than in 
the shallow-water equations. 

Figures l-[6] show examples of Solar System giant planet simulations from our study. 
The global planetary parameter values adopted for the simulations are listed in Tabled! The 
values for HD 209458 b are also included in the Table for comparison. In the simulations, 
small finite eccentricities are ignored. Since it is known from observations that the emitted 
flux is fairly uniform over the entire surface of all Solar System giant planets, we do not 
include any differential thermal forcing in our simulations of these planets. The qualitative 
success of these calculations in reproducing the global circulation pattern of Solar System 
giant planets gives confidence in our approach to modeling these planets and unsynchronized 
extrasolar giant planets with small eccentricities. 

Figures [1] and [2] show a typical evolution from one of our Jupiter simulations (Model SI 
in Table [2]). Contour maps of the potential vorticity (g) field are presented. Recall that q 
is a tracer of the fluid, hence a proper variable to plot. Positive values are contoured with 
full lines, and negative values are contoured with dashed lines. In Figure (H the maps are 
in orthographic projection centered on (A, 0) = (270, 0) and gridded at 30-degree intervals. 
Note that with this projection the polar latitudes occupy a smaller area of the disk than 
the low latitudes. The resolution of the simulation is T170 — i.e., iV = 170 in Eq. (|T2|) — and 
corresponds to a 512x256 longitude-latitude grid over the globe. This resolution is minimally 
comparable to a 1024x512 grid in a finite difference calculation because of greater intrinsic 
accuracy of the pseudospectral method. The duration of this run is 2000 planetary rotations 
(i.e., r = 2000). We have also performed shorter duration runs at higher resolutions (T213 
and T341) to verify convergence of the results and find that there is no qualitative difference 
in the evolution at the higher resolutions. However, a minimum of T63 (192x96) resolution 
is required to resolve the turbulence (i.e. growing vortices). Time in unit of the planetary 
rotation periods is indicated in the upper left corner of each frame. The contour levels, 40 
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in all, are identical in each frame. In Figure [2], the maps are in polar stereographic views 
centered on the north pole. In this projection also, the polar region occupies a smaller area 
of the disk than the equatorial region. East is in the counter-clockwise direction, and A = 
is at the 3 o'clock position. 

At r = 0, the atmosphere is stirred randomly and released with a background correspond- 
ing to the planetary vorticity /. By r = 84, the formation of zonally elongated structures 
is clearly evident at low latitudes. These structures are formed by growing eddies, which 
radiate Rossby (planetary) waves after reaching the size of Lp (see §2). Rossby waves are 
the undulations of q about a latitude circle. Note that the wave propagation is strongest at 
the equator since the linear dispersion relation from Eq. (|TT|) for the Rossby waves gives a 
longitudinal phase speed, 

^ ~n{n + l)/^l+l/LnV ^^^^ 

Recall that (3 = R~^df /dcf) = 2Q cos (p/Rp and is a maximum at the equator. The length, 
Ltis = {gHeY^"^ / f , is the "equivalent depth" Rossby deformation ra dius, wher e is the 



equivalent depth — an ir nportant parameter in atmos pheric dynamics (jGilllll982l ). including 



in the theory of tides (1 Chapman and LindzenI llQTOl ). In the equivalent-barotropic model 



He = Hpo/Ao, a weighted pressure scale height evaluated at the lower bounding surface. Note 
also that is largest at the equator and smallest at the pole, giving a small interaction 
length for the vortices in the polar region (see r = 423 frame in Figure [2|, for example). The 
small interaction length significantly reduces motion and mergers in that region. From our 
simulations, we expect the polar region at the cloudtop level of Jupiter to contain spots, 
rather than bands. And, in fact, bands are not observed above about 60 degrees latitude on 
Jupiter. 

Maps of q can be effectively used to study atmospheric circulation. Such ma ps constitute 



a standard tool for studying the circulation on the Earth (IHoskins et al.lll985l ). They have 
frequently provided insights to important processes which are otherwise difficult to obtain. 
Since q = g(v), under certain "balanced" conditions (e.g., pressure gradient force balanced 
by the Coriolis fore . q can be inverted to obtain the wind field, v = v(g). That is, the 
wind information is contained in q. In spectral space, q can be seen to have a broad distribu- 
tion. This demands a high resolution for accurate representation, both observationally and 
numerically. Loss of q due to under-resolving, therefore, causes loss of pertinent information. 
Indeed, the ability to retain q is one strong justification for using the equivalent-barotropic 



^^This condition is known as geostrophic balance ( Holtonlll992[ ) and applies to Solar System giant planets, 
as well as the Earth, to a good degree away from the equator. 
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model and also an important advantage of our model compared with models using the finite 
difference method. 

Finely resolved q features can be used to identify the presence of atmospheric flow 
stru ctures (e.g. , vortices and jets), as well as flows around the structures, to exquisite detail 



(cf., ICho et al.l (120011 ). and references therein). For identifying jets, consider a simple 2-D 
zonal flow: u = u{(j)) and f = with u, f G 3?^. If incompressible (non-divergent) in 2-D, Ti. 
is constant. Then, Dq/Dt = with q = (C + /), and TC plays no dynamical role. Using the 
definition of ( (= k-Vxv), we have ( = —R~^du/d(j). Note that if i?o <S 1 and ( ^ U/L, then 
C -C /. Therefore, for rapidly-rotating planets (such as Jupiter), the Coriolis parameter / 
dominates over the relative vorticity (, especially at high latitudes. In this case, q is mostly 
positive (negative) in the northern (southern) hemisphere, as can readily be seen in Figured! 
Recall that / = 2Q sin and hence changes sign when crossing the equator. The gradient of 
q in the meridional direction is then 

where U = ucos(f) and /i = sin0. On the plane tangent to the planetary surface, this is just 



- = /?-—, (21) 
dy dy"^ ' 



where [3 = df /dy and y is the northward direction. According to Eqs. fl20l) and fl2T|) . jets 
may be identified by their curvatures. For example, in Eq. fl2Tl) . an eastward jet corresponds 
to d'^ujdy^ < while a westward jet corresponds to d'^u/dy'^ > 0. Consequently, for a given 
value of f3 (which is always positive), a large (positive) value of q gradient indicates an east- 
ward jet in the fiow while a small (negative) gradient of q indicates a westward jet. Thus, an 
atmosphere containing alternating eastward and westward jets show a corresponding alter- 
nating tightening and relaxing of g-contours in the meridional direction. The correspondence 
can be seen in Figures [2] and [3l The absolute value of dq/dy (or dq/dcj)) is a direct measure 
of the sharpness of the jet, eastward or westward. 

Figure [3] depicts the quasi- steady jets from the Jupiter simulation presented in Figures[T] 
and [2] (Model SI). It shows a zonally-averaged (eastward) wind profile as a function of 
latitude, [u], where 

K0)] = 7^ [\{<l>,X)dX. (22) 
27r Jo 

The jets are the peaks in the wind profile. The time of the simulation is r = 300, well 
after the profile has emerged from the initial stirring. The profile is steady and nearly 
identical even at r = 2000. In the figure, there are about half dozen zonal jets in each 
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hemisphere with amplitude ~25 m s~^. There is also a very strong equatorial jet of amplitude 
~ —125 m s~^ Qualitatively, the number, widths, and absolute magnitudes of jets are all 
similar to those observe d on Jupiter. As in the previous shallow-water studies of Jupiter 



( ICho and Polvanilll996bl ). however, the sign of the equatorial jet is opposite to that observed 
on Jupiter, as well as on Saturn (not shown). The retrograde (westward) equatorial jet in 
the simulation is due to the strong presence of Rossby waves at the equator, as discussed 
above. Clearly, an additional driving mechanism is also required in the equivalent-barotropic 
equations model to capture the proper sign of the equatorial jet, as in the shallow-water 
equations model. This is not unexpected given the formal similarity between the equivalent- 
barotropic equations and the shallow-water equations. Later, in §5, we return to the issue of 
prograde (eastward) equatorial jet — currently an unsolved problem in planetary atmospheric 
dynamics. For now, we simply note that the model does capture the proper sign of equatorial 
jets on Uranus and Neptune, as we shall show shortly. Evidently, some mechanism not 
inclu ded in the simp l est ve rsion of the equivalent-barotropic equations is required for Jupiter 



(cf., iHeimpel et al.l (120051 )) and Saturn but not for Uranus and Neptune. We also note 
that simulated bands are more wavy and contain fewer spots than on the observed planets. 
However, these simulated features appear to match the observed ones better with increased 
resolution in the simulations. 

In Figure HI we present a direct comparison of the equivalent-barotropic and shallow- 
water models. We do this to demonstrate the validity of our model. When k = 1, the adiabatic 
equivalent-barotropic equations formally reduce to the shallow-water equations, with 
set to unity0 In this case, HpQ becomes the layer thickness in the shallow-water model 
(Eq. (Q). Note that, while formally equivalent, the "shallow-water limit" of equivalent- 
barotropic equations does not carry a physical interpretation. This is because physically 
K = 1 corresponds to c^, = and 7 — oo, where is the specific heat at constant volume 
and 7 = Cp/c^ (recall that k = R/Cp). In the figure, two runs are presented: on the left 
(a) is a frame from the equivalent-barotropic calculation and on the right (b) is a frame 
from the shallow- water calculation (Models Jl and J2 in Table [2], respectively). In the two 
runs, all parameters are identical, except that the equivalent-barotropic equations model 
calculation on the left has k = Ao = 1. The resolution of the runs is T106. The projections 
of the maps are polar stereographic, as in Figure [21 The time of the frame is at r = 179, 
long after the quasi-steady state has been reached. In these two runs, the initial spectrum is 
peaked at a lower wavenumber (larger scale) than in the run presented in Figures [1] and [2l 
As a consequence, the steady state is reached earlier in these runs, compared to the run in 



^^We remind the reader that vortex columns cannot tilt in the shallow-water model. So, a unit ^0 is 
physically consistent. 
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Figures [T] and [2l However, because both runs start with peaks at scales much smaller than 
the Rhines scale Lp, there is no qualitative difference in the final flow configuration: only 
the time when the equilibrium state is reached is different. Significantly, the results from 
the two models, (a) and (b), are practically identical0 This verifies that flow features are 
truly robust (independent of the randomness of stirring, for example) and are reproducible 
by the new model. 

Figure shows three successive equatorial orthogonal (a) and polar stereographic (b) 
views of a typical Neptune simulation from our model (Model S4 in Table [2]). In the figure, 
three time frames are shown to illustrate i) the flow in the initial state (r = 0), ii) during 
transition to the quasi-steady equilibrium state (r = 33), and Hi) the equilibrium state itself 
(r=146). As in Figures [T] and [2], iso-g contours are plotted with positive (negatives) values 
in full (dashed) lines. The initial spectral distribution is the same as in that of Figures [T] 
and [21 except U is larger (see Tabled]). In this simulation, the structures dominated by 
scales around the wavenumber n ~ 14 at r = are more easily recognized than in the Jupiter 
case because of the larger U and smaller / for Neptune. More importantly, because of the 
larger U and smaller /, the flow quickly evolves into a final configuration dominated by two 
strong circumpolar vortices (one centered on each pole) and a broad jet, associated with a 
well- homogenized q region at low latitudes. In fact, only a few g-contours are seen all the 
way up to the high latitudes, until the boundaries of the polar vortices are reached. Such a 
q = low-latitude region was previously assume d in modeling the ch aotic, wobbling motion of 



Neptune's Great Dark Spot with great success (iPolvani et al.lll99Cll ). The strong polar vortex 
results from vigorous mergers, in which surrounding vortices are continuously "soaked up" 
before they are sheared away or dissipated. Dynamically, the strong interaction is expected 
given the large value of L-ji/Rpr^ 1/6 for Neptune, in contrast with ~l/30 for Jupiter (see 
Table [2]). Note that the circumpolar vortices are cyclonic. Under geostrophic balance, the 
cyclonic vortex is associated with a cooler region than its surroundings. 

In the case of Neptune as well, the steady-state configuration (already well in place 
after r ~ 30) is stable over timescales of several thousand r. Significantly, as will be shown 
shortly, Neptune's jet profile is closer to that we obtain for close-in extrasolar giant planets 
than Jupiter, as reported in ICho et al. J2003h . This marked difference mainly arises from 



the lower VL values of close-in extrasolar giant planets compared to that of Jupiter. In 
fact, assuming synchronization, close-in extrasolar giant planets rotate only moderately fast 
compared to Solar System giant planets. Additional differences between close-in extrasolar 
giant planets and Neptune also result from the intense irradiation experienced by the close-in 



^^They of course cannot be exactly identical because of the initial random seeding and minor numerical 
differences. 
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extrasolar giant planets, as we shall see below (§5). 

But, first we quickly summarize our equivalent-barotropic equations results for the Solar 
System giant planets. Figure M shows the zonal wind profiles at late times obtained in our 
simulations of all four Solar System giant planets (Models S1-S4 in Table[2]): Jupiter (a), Sat- 
urn (b), Uranus (c), and Neptune (d). Qualitatively, these profiles match well the observed 
profiles in terms of the overall number of bands and strengths of jets on each planet (see Cho 
& Polvani 1996b for observed profiles). Generally, the four planets organize themselves into 
two discernible groups in terms of their general circulation patterns at the cloudtop levels. 
Jupiter and Saturn share similar fiow patterns, with their multiple jets and narrow bands; 
and, Uranus and Neptune share similar fiow patterns, with their few jets and broad bands. 

In the current study, we have performed over 300 simulations in order to thoroughly 
explore the numerical and physical parameter-space spanned by the equivalent-barotropic 
equations system. Based on this extensive exploration, it appears that capturing the pro- 
grade rotation on Jupiter and Saturn requires physics beyond that included in the current, 
adiabatic equivalent-barotropic equations model. We give one example of a possible ex- 
tension in §5. It is crucial to note that the lack of correspondence is not simply an issue 
of including the third dimension. One well-known me chanism in 3-D situatio n which does 



produce prograde rotation in simulations of the Earth (ISuarez and Duffylll992l ) is a zonally- 
asymmetric (m = l and m = 2), large-amplitude forcing that is symmetric in the meridional 
direction. However, this mechanism does not apply to Jupiter (since it is not thought to 
be zonally asymmetrically heated) and cannot explain Saturn (since it is not meridionally 
symmetrically heated). In §5, we discuss the consequences of this forcing further for close-in 
extrasolar giant planets, where it may find some justification. 

In summary, the results from our study illustrate several generic features of giant planet 
circulations modeled with the equivalent-barotropic equations. First, there is no qualita- 
tive difference between model results for equivalent-barotropic equations and shallow-water 
equations. Second, with the knowledge of a few global physical parameters, we are in a po- 
sition to capture the gross properties of unsynchronized extrasolar giant planet atmospheric 
fiows (several examples of possible unsynchronized extrasolar giant planet circulations are 
presented in the ensuing sections). Third, as in the previous shallow- water and equivalent- 
barotropic studies, the scales Lp and L-ji emerge as critical dynamical parameters for effecting 
realistic "predictions" of general circulations. Last, a knowledge of the strength and distri- 
bution of diabatic forcing (i.e., not just incident but also barochnic, convective, latent heat, 
etc.) appears minimally needed to properly model the direction of the equatorial jet. 
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4.2. A Warm Jupiter 

Figure [7] illustrates our simulation of a "warm" Jupiter (Model Ul in Table [2]). This 
may serve as an example of a Jupiter-like extrasolar giant planet which has migrated in 
some distance toward its host star but is not yet synchronized. At the same time, since the 
temperature increases with depth below the main cloud deck on Jupiter, it may heuristically 
serve as a representative flow at a deeper level in the Solar System planet. In this simulation, 
all the parameters are identical to that shown in Figures [1] and [2] (Model SI), except Ti in 
this simulation is such that the global average temperature To = 300 K (roughly a "Jupiter 
at 1 AU"). Both stereographic (a) and orthographic (b) views are shown at r = 199, well 
after the formation of steady bands and jets. This figure should be compared with the latter 
time frames of Figures [1] and [2] (Model SI). 

In the early part of the evolution, the flow is very similar to the early evolution of Model 
SI. However, small differences begin to appear as the evolution proceeds: the banding at 
low to mid latitudes is weaker than that in Model SI and even shows signs of slight erosion 
at times. The effect of the erosion can be seen in Figure [TJ where the positive gradient of 
q at low latitudes is weaker compared with those in Figures [1] and [21 The result is that 
bands generally appear more diffuse and there are one or two fewer jets in each hemisphere. 
The amplitude of the equatoria l jet is reduced by more than 50% as well. This behavior is 



consistent with the findings in ICho and Polvanil (jl996al ) for the shallow-water case and is 



due to the increase in leading to a larger (Rossby deformation radius). The larger 
L-ji acts as increased "stiffness" , or reduced horizontal compressibility; this leads to stronger 
interaction between the structures in the fluid and the bands themselves undergo several 
"mergers", reminiscent of the vortices. Formally, the shallow-water equations reduce to 2-D, 
incompressible equations in the limit of L-j^-^ oo. Note that in the figure no mid-latitude 
anticyclonic spots survive and only a few vortices are present in the polar region, as a further 
consequence of the larger L-jz. 

Given the behavior illustrated in Figure [TJ we may draw several general conclusions. 
First, a much closer-in Jupiter may have a slightly different visual appearance than the 
present Jupiter at 5 AU. It will possess broad, diffuse bands with very few, possibly no, 
spots at the main cloud deck. Second, except for the equatorial jet, the overall jet structure 
is not changed very much, attesting to the strong control rotation has on the flow. The ratio, 
L-ji/Rp, is still <^ 1/3 in this simulation and a complete breakdown of the basic jet structure 
is not expected: an unrealistically large Tq ^ 1.3 x 10^ K is necessary for the breakdown. 
Finally, given that the situation can also be considered as crudely modeling a deeper level 
(one pressure scale height lower) in the current Jupiter, the flow is fairly vertically coherent 
(i.e., barotropic). There is not much difference in the flow compared with that at the "higher 
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level" (i.e., flow with Tq = 130 K), even in the absence of the explicit vertical coupling present 
in the baroclinic models. 



4.3. Unsynchronized HD 209458 b 

Having established the singular importance of the rotation rate Q, by considering 
Jupiters with two different equilibrium temperatures (or orbital radii), we now consider 
the case of a "Jupiter" when its Q is signiflcantly reduced (e.g., by tidal interaction) but still 
not locked in 1:1 spin-orbit resonance. Such a planet would have a small semi-major axis, 
a ^ 0.2 AU. This is an important case to study since it allows a clear delineation of the con- 
sequences of tidal locking on the atmospheric circulation. It i s also relevant for the rec ently 



proposed scenario for HD 209458 b caught in a Cassini state (jWinn and HolmanI 120051 ). As 
such. Figure 6 presents one of the key results of this study: atmospheric circulation on a 
planet with parameters representative of HD 209458 b (see Table [T]) but without the explicit 
asymmetric forcing expected in the case of tidal synchronization. 

In Figures [8] and [9|, the flow evolution from Model U2 (see Table [2]) is shown at the early 
stage (a) and at the equilibrated stage (b), respectively. Once the latter stage is reached, 
the flow state does not change in its essential character: the circulation is dominated by 
a single, coherent vortex at the pole and a sharp jet at its periphery, reminiscent of those 
on Neptune/Uranus rather than Jupiter /Saturn (compare with Figures [2] and [5]) . As in the 
simulations of Solar System giant planets, the initialization here is a random stirring, which 
is nonlinearly balanced and consistent with a specifled Tq and U. In this simulation, Tq = 
1450 K and U = 400 m s~^, plausible values for HD 209458 b. The actual values, particularly 
U, are currently unknown. Note that these values give L-ji/Rp ~ 1 and Lp/ Rp ~ 1, as shown 
in Table H 

Although clearly more similar to the circulation of Neptune/Uranus, there are some 
distinguishable differences. The polar vortex in Figure [HI is much larger than in the Neptune 
case. Also, the vortex is not centered on the pole, as it is on Neptune (see r = 146 frame 
in Figure [5]). Instead, it can be seen revolving around the pole on a timescale of about 2.7 
planetary rotation periods in this run. In general, the timescale of the revolution around 
the pole depends on the strength and position of the polar jet (related to f7), since the 
polar vortex essentially moves with the jet (background zonal flow). In our extensive survey, 
we have found that the timescale of the revolution around the pole range from a few days 
to several tens of days. In addition, comparing the r = 146 frame in Figure [5] with any 
of the frames in Figure [9], one can identify a marked reduction of Rossby wave amplitudes 
(undulations of q lines about a latitude circle) at low-latitudes in the extrasolar planet 
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case — although wave breaking (lateral overturning of a large-amplitude wave) occurs on 
both planets. Associated with this feature is the reduced mixing and homogenization at low- 
latitudes of this planet, compared with Neptune. All of these behaviors can be understood in 
terms of the larger Ln and smaller f5 on the unsynchronized close-in extrasolar giant planet. 
Larger L-ji means more robust mergers and a more dynamic final vortex. Smaller f3 entails 
a weaker shearing environment for vortices at low latitudes and smaller wave amplitudes 
there. 

An off-the-pole polar vortex, as in Figure [9l has several significant consequences for 
mixing of heat and vertical motion at high latitudes. If the vortex is cyclonic, as in the 
figure, the air in its interior is cooler than it's surroundings. Hence, a translating vortex 
carries with it a large mass of cold air around the pole. Moreover, as the vortex translates, 
it pulls up warm air from the lower latitudes ahead of it while pulling down cold air from 
higher latitudes behind it. This induces a poleward flux of heat and meridional mixing. The 
effect is accentuated by its meandering motion about a latitude circle. At the same time, 
the vortex cools in the net because it is generally warmer than the radiative equilibrium 
temperature (recall that very little stellar irradiation reaches the high latitudes). Hence, 
by the first law of thermodynamics, the potential temperature of the air inside the vortex 
decreases. But, since positive stability of a stratified region implies 89/ dz > 0, where 9{z) is 
the basic potential temperature, the cooling air must sink across ^-surfaces to a lower valued 
surface; that is, there is a slow, down-welling motion inside the vortex (not modeled in this 
work). In addition, planetary wa yes further i ncrease the downward drift by a hysteretic 



mechanism in diabatic conditions (ISalbylll996l ). Later, in §5, we will show that planetary 



waves are induced by thermal forcing on extrasolar planets. 

Figure [TOk shows the quasi-steady jet profile from the simulation in Figures [8] and [9] 
(Model U2). In this simulation, in which U = 400 m s~^ has been adopted, three broad 
jets form. The low number of jets is also a robust feature of our HD 209458 b simulations 
with day-night forcing (§5.1). While the magnitude and sign of individual jets vary from run 
to run, the number and root mean square speed of the jets mainly depend on U, in rough 
agreement with L^. Occasionally, when U ^ 400 m s~^, a two-jet profile forms. The adopted 
f7 value is characteristic of Saturn, Uranus, and Neptune (see Table [T]). On extrasolar giant 
planets, U can plausibly vary from ~100 m s~^ to ~1000 m s~^. A value several times 
greater than ~1000 m s~^ is not likely to be realistic for large-scale flows on extrasolar giant 
planets since the sound speed Cg ^ 3000 m s^^ (~2700 m s^^ on HD2O9458b)0 The two-jet 
configuration also forms in the special case when U = 0, as we will show in §5.3. At the low 



^This need not be so for small-scale waves, which may propagate with speeds close to or at Cg. 
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end of the range, when U = 100 m s ^, a four or five jet profile can result, consistent with 
L^. This case is illustrated in Figure [TOb. 

At this point we wish to remark on an issue pertaining to jets which requires clarifica- 
tion. There appears to be some confusion in the current literature concerning the universality 
of retrograde equatorial jets in one-layer turbulence calculations. We emphasize here that a 
retrograde equatorial jet is not a necessary outcome of all shallow-layer tur b ulence calcula- 
tions. That outcome requires a small L-n, as pointed out in lCho and Polvanil (Il996al ). In the 
case of Solar System giant planets, L-ji is indeed small for all the planets, and the simulated 
equatorial jet is always retrograde in the absence of forcing. On extrasolar giant planets, 
L-ji need not be small. In fact, L-ji is generally not small for close-in extrasolar giant planets 
because they are much hotter and less-rapidly rotating than Solar System giant planets. 
Accordingly, the equatorial jet in Figure fTOb is actually prograde. We have performed many 
simulations of HD 209458 b under a variety of conditions and find that what is robust is the 
adherence of the general flow pattern to Lp, not the direction of any particular zonal jet. 

As already discussed, another dynamically relevant factor that is currently unknown 
is the radiative equilibrium temperature distribution on extrasolar giant planets. The dis- 
tribution is related to the stellar irradiation but also depends in a complicated way on the 
composition, cloud and aerosol distribution, and the dynamics itself through the role it plays 
on the distribution of the radiatively active components. Figure [TT] shows the simulation pre- 
sented in Figure 6 (Model U2) at two different equilibrium temperatures: (a) 1800 K and 
(b) 800 K (Models U3 and U4, respectively). The range of the two temperatures is broad and 
is meant to roughly bracket Tq resulting from a variety of physical situations (e.g., distance, 
stellar flux, and albedo). With the same flux and albedo assumed in Model U2, the higher 
To can also be applied to the case of an HD 209458 b-like planet which is ~0.005 AU closer 
to the host star than HD 209458 b at present (a = 0.045 AU). Similarly, the lower Tq can 
be applied to an HD 209458 b-like planet that is ~0.15 AU further away. As can be seen in 
Figure [m the higher Tq leads to a more smoothly varying low- to mid- latitude q distribution 
while the lower Tq leads to much more planetary (Rossby) wave activity, along with a more 
centered polar vortex. In both cases, the basic three-jet profile is unchanged. However, 
interestingly, the higher Tq value leads to di splaced a polar vortex, which may potentially be 
of observational relevance (ICho et al.ll2003l ). 
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5. Synchronized Extrasolar Giant Planets 

5.1. Effects of Zonally Asymmetric Heating 

From an atmospheric hydrodynamics point of view, HD 209458 b is particularly interest- 
ing in that both its radius and mass (hence its surface gravity) are accurately known and it 
has been the subject of most observational studies to date. The radius and mass information 
is now available for more than a dozen extrasolar giant planets and this number is expected 
to grow in the future. As already noted, even if synchronization is invoked to pin down Q 
and a "cloudless atmosphere" is considered to leave aside the complexities of radiative issues 
(clouds, albedo, vertical distribution of species, etc.), there is still no information available to 
constrain the amplitude of thermal forcing rj and the global wind speed U. These parameters 
crucially control the flow dynamics through local changes of L-ji and global changes oi Lp. 
Eddy/wave amplitudes and propagations are affected as well. Experience with Solar System 
planets tells us that circulation and radiative transfer mutually interact in such a way as to 
modify the temperature structure away from the simple radiative and radiative-convective 
equilibrium. In this section, having considered a "HD 209458 b" in an unsynchronized state, 
we now isolate the effects of a specifled value of r] on the flow. In §5.21 we present our flndings 
of the effects of varying f7 on the flow and the horizontal temperature structure. We wish to 
focus on the robustness of plausible circulation patterns formed under turbulent conditions. 

As a flrst study of the forcing effect, we assume a simple day-night equilibrium tem- 
perature distribution, z'^ oc (cos A cos 0)^*, where x is chosen from the set of {1/4, 1, 2} for 
different runs. Note that a full 1-D radiative transfer calculation of HD 209458 b atmosphere 
applied over the planet suggests that this is a good range for the planet, with x= 1/4 being 
more appropriate for the upper levels and x=l being more appropriate for the lower levels 
in the stable layer of a close-in extrasolar giant planet. There is no qualitative difference in 
our results due to a different choice of x value adopted, however. The main result of our 
study of the forcing effect is as follows: in adiabatic, equivalent-barotropic, turbulent sim- 
ulations, the global circulation is not qualitatively affected by the applied thermal forcing. 
That is, the general characteristics of global circulation (i.e., number of bands, strength of 
zonal winds, presence/absence of dominant polar vortices) are still well characterized by L-ji 
and Lp (or equivalently Ro and 5„), irrespective of the value of r] adopted. This may appear 
surprising at flrst, but one of the key points of the present work is that the basic flow pattern 
is set by the dynamical parameters, Lp, L-ji, and Rp. Unless very strong diabatic forcing 
places the extrasolar planet into a special dynamical regime (see §5.3), forcing which merely 
deflects the isentropic surfaces (as opposed to causing large mass to cross the surfaces, as 
in overturning surfaces) cannot in itself produce major changes in the flow at this level of 
modeling. It can, of course, affect the temperature structure, as will be shown shortly. 
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The robustness of the basic flow is illustrated in Figure [121 The figure presents four 
T106 runs (Models HI, H2, H3, and H4 in Table [2]), which are identical in all respects except 
for the forcing perturbation amplitude, zq. In the figure, zq spans from to 20% of the 
background temperature (i.e., rj = to 0.2), corresponding to to ~ 300 kH The latter is 
not a small deflection. We have also run cases in which the fractional deflections are even 
larger (77 = 0.4), but the behavior is basically a more violent version of the 20% case. Frames 
at the same time (r = 95) from simulations with rj values of (a), .02 (b), .04 (c), and .20 (d) 
are shown. As seen from the figure, the forcing primarily amplifies the mid-latitude jet. As 77 
increases, low latitude q is homogenized, signifying strengthening of the retrograde equatorial 
jet (a-c). Concurrently, the g-gradient at mid-latitude, hence the jet there, becomes sharper. 
The wave breaking activity in the mid-latitudes also diminishes, until the forcing is no longer 
small (a-c). In the large-amplitude case (d), the jet becomes strongly zonally asymmetric, 
becoming unstable, along with increased wave breaking activity in the mid-latitudes. When 
dq/dfi changes sign, due to the sharpening of the jet, t he jet meets the necessary condition 



for barotropic instability and the unstable flow ensues (jCharney and Sternlll962l ). 



Figure [13] shows the effects of increased forcing amplitude 77 on the temperature distri- 
bution. To = ro(A, 0), at dynamical equilibrium. The same runs in Figure ITU (Models HI, 
H2, H3, and H4) are shown with 77: (a) 0., (b) .02, (c) .04, and (d) .20. As can be seen, for a 
fixed If, Tq is strongly dependent on rj. With no or low forcing, Tq essentially tracks the flow, 
the polar vortex (a). As the forcing amplitude increases, however, competition between the 
flow-induced moving temperature pattern and the forcing-induced fixed temperature pat- 
tern develops: temperature anomalies on the whole either rotate quasi-steadily about the 
pole (b) or oscillate aperiodically about a point near the substellar point (c). When the 
amplitude is high, the forcing pattern dominates entirely and Tq is essentially the imposed 
day-night distribution fixed on the substellar point (d). We emphasize here that the global 
circulation is not modified much, at least in the adiabatic case. Similar jet profiles and polar 
flow structures are present, even under strong forcing (compare Figure [12] with Figure [T^ : 
an asymmetric^ cyclone/anticyclone pair at high latitudes is always present (see Figures [T2l 
and [15]), but it may not be discernible in temperature maps if swamped by the applied ther- 
mal forcing. We note that when the same set of runs is repeated with larger U, the same 
trends in temperature maps are exhibited at comparatively larger rj values, as expected from 
stronger motion-induced temperature patterns. Hence, thermal variability depends on U, 
as well will be further demonstrated in the next subsection. This is one of our key 

results. 



^^Recall that the average background is represented by average thickness (mass) of the layer. 
^"^The cyclone is much stronger than the anticyclone. 



- 36 - 



In summary, clearly, the entire spectrum of temperature distribution behavior is possible — 
from forcing being completely overwhelmed by the flow to forcing completely overwhelming 
the flow. Not unexpectedly, the critical t] value at which Tq makes a transition to merely 
expressing the specified forcing depends on the given U value. The fact that the flow is fairly 
similar in all these cases indicates that the coupling between the flow and applied heating 
is relatively weak. Diabatic forcing appears to be required for a stronger coupling. The 
overall implication here is that, one cannot infer U and t] independently of each other from 
a temperature map, at least in adiabatic conditions. 



5.2. Global Mean Kinetic Energy Variation 

Figure [14] shows the dependence of the global circulation on U. The figure presents 
four T106 runs (Models HI, H5, H6, and H7 in Table [2]) that are identical in all respects, 
except for U. In these runs. If spans from 100 to 1000 m s~^ from run to run. As noted 
already, the latter is not a small value: globally-speaking, the circulation is no longer in the 
quasi-geostrophic regime, when U ~1000 m s~^, since Ro^l. Locally, of course, Ro may be 
<^ 1, since the velocity varies over the planet. In the figure, simulations with U values of 100 
m s~^ (a), 200 m s~^ (b), 400 m s~^ (c), and 1000 m s~^ (d) are presented at the same time 
frame (r = 95). We have also run cases in which U is even larger (2000 m s~^), but in this 
case the behavior is similar to that of (d), at early times. At later times, the flow field blows 
up, due to the large amount of surface gravity waves that naturally ariseQ It can be seen 
clearly in Figure [H] that increasing U increases the size and strength of the polar vortex. In 
fact, the vortex is large enough in (d) to nearly cap the planet poleward of 60 degrees if it 
were centered at the pole. This results from the greater root mean square kinetic energy, 
|f7^, of the initial eddies. 

It may appear to the reader that, in Figure [HI planetary (Rossby) wave breaking activity 
abruptly terminates at high U. However, in actuality, it increases. This is demonstrated 
in Figure [TBk-c. in which the field presented in Figure UM i (reproduced in Figure fTEk) is 
successively contoured at smaller range and intervals (Figure [T5b-c): the number of contours 
is the same. The vortex is "chopped off" to reveal the surrounding flow. Clearly there is a 
large amount of structure outside the dominant vortex. In fact, the large amplitude of the 
planetary (Rossby) waves and their breaking has produced an anticyclone (Figure [TBb). This 
anticyclonic vortex corresponds to a warm region compared to the interstitial region between 



^^The blowing up can be prevented, to a certain extent, if an unrealistically large amount of artificial 
dissipation is applied or, equivalently, if the resolution is very low. 
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the pair, giving rise to an associated slowly-rotating thermal dipole. This is remarkably 
similar to what happens on the Earth: in the Earth's stratosphere, upwardly-propagating 
planetary waves cause the stratospheric polar vortex to be displaced from the pole and induce 
the formation of an anticyclone, conserving q. The interstitial region itsel f is well mixed due 



to pronounced wave breaking, as on the Earth (cf., [Andrews et al.l (Il987l )). 



Figure [TSti-e shows the effect of increasing IJ on the zonal winds. The zonal wind 
profile is plotted for two simulations presented in Figure [HJ (a) and (d) — Models HI and 
H7, respectively. The figure clearly shows that, for fixed VL and Rp, the number and width of 
the jets depend on U in good agreement with the Lp estimate. Compare also this figure with 
Figure [T0|, in which asymmetric heating was not applied. From the comparison, one may 
draw two general conclusions: 1) forcing does not appear to have a distinguishable effect on 
the general structure of the zonal jets, their number and width; and 2) given the broad range 
of U values in the simulations of the two figures, the low number of zonal jets is a robust 
feature on HD 209458 b-like planets. 

The tJ values can be shown to be intimately tied with the resulting temperature distri- 
bution, as shown in Figure [TBI The same runs as in Figure [T3] (Models HI, H5, H6, and H7) 
are shown with U values: (a) 100 m s~^, (b) 200 m s~^, (c) 400 m s~^, and (d) 1000 m s^^. 
As can be seen, for a fixed forcing amplitude ?7, Tq is strongly dependent on U . This is 
very similar to the behavior illustrated in Figure [131 With small U , Tq essentially expresses 
the day-night difference (a); that is, the flow is not energetic enough to advect tempera- 
ture in significant amounts. As U increases, however, competition between the flow-induced 
temperature pattern and the forcing-induced temperature pattern develops: temperature 
anomalies on the whole either oscillate aperiodically about a point near the substellar point 
(b) or rotate quasi-steadily about the pole (c). When U is high, the flow dominates and Tq 
essentially reflects the distribution dictated by the dominant vortices (d). Again, in each 
case, there is an asymmetric cyclone/anticyclone pair at high latitudes, independent of U . 
When the same set of runs is repeated with larger rj values, the same trend is exhibited, with 
the transition to flow dominance occurring at larger U value. The only time a qualitative 
departure is obtained is in the special case when the simulation is started at a resting state 
and driven wholly by the applied day-night forcing (without any stirrings), the situation to 
which we now turn. 



5.2.1. Quiescent Initial State: A Special Case 



At the current level of modeling, be it with the equivalent-barotropic equations or the 
full primitive equations, one of the ingredients necessary for accurate modeling and good 
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understanding of extrasolar giant planet atmospheric circulation is the physics of baroclinic 
processes. To be more precise, what is important is the inclusion of the effects of the 
processes, not so much the actual processes themselves. This is the idea of parameterization. 
In our usual model setup, the atmosphere is stirred in the beginning. The stirring is designed 
to crudely represent a generic form of baroclinic or barotropic processes: eddies and waves 
generated by small-scale convection or baroclinic instability, for example. Both of these 
generating mechanisms are expected to exist on extrasolar giant planets. The thermal forcing 
should be another source of strong eddies and waves, especially on close-in extrasolar giant 
planets. Once created, the eddies and waves tra nsport and/or propagate crucial dynamical 
quantities, such as heat and wave-activity densitjo ( Andrews et al.lll987l ). For example, the 
action of waves greatly affect both the source and the target region (where the waves are 
dissipated) via interaction with the background flow (e.g., jets) present in the target region. 
This is one reason why we have emphasized resolution and laid particular stress on the jet 
profiles, and flow structures in general, in this study. 



Recently, IShowman and GuillotI (120021 ) and lCooper and Showman! (120051 ) have performed 
baroclinic simulations of close-in extrasolar giant planet atmospheres using the primitive 
equations. In these simulations, the atmosphere initially at rest is driven by day-night forcing 
applied in the energy equation. As in some simulations of "stripped-down" Earth, which em- 
ploy the primitive equations and similar forcing (i.e., one that is low-order and meridionally- 
symmetric, but zonally-asymmetric), they obtain close-in extrasolar giant planet atmo- 
spheres which are superrotating. That is, the zonal wind at the equator is prograde: 
[m(0 = O)] > 0. In Earth simulations, the emergence of superrotation depends nonlinearly 
on the amplitude of the forcing, sometimes ma king an abrupt transition from subrotating to 
superrotating state at some critical amplitude (ISuarez and Duffylll992l ). At criticality, a bi- 
furcation is observed, the flow exhibiting sometimes superrotation and sometimes not, for an 
identical set of parameters. In some cases, the superrotation appears to be in a stable equihb- 
rium state, as the atmosphere remains superrotating even when the forcing is removed. Ex- 
actly how this peculiar behavior arises is currently not well-understood. Feedbacks between 
eddy momentum, angular momentum, Hadley cell, bottom friction, and vertical resolution 
appear to be at play. Certainly, vertical angular momentum exchanges between atmospheric 
layers and horizontal exchanges within a layer is required for superrotation. 

In simulations of the Earth, superrotat ion has also bee n achieved, and investigated, 
through direct angular momentum forcing (jSaravananI Il993l ) . As with many phenomena 
admitted by the full primitive equations, superrotation of planetary atmospheres can also 



^^For homogeneous fluids, wave-activity density is [^(^']/[dq/dij], where (' is the vorticity-eddy. It is 
sometimes also called (minus) pseudomomentum in the literature. 
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be s tudied using the equ ivalent-barotropic equations in a similar way. In this work (see 
also lShell and Heidi (120041 )). an equatorial torque is applied to represent what is expected to 
happen to and within an individual layer in a (multi-layer) baroclinic atmosphere. This is 
done by adding a forcing term, 

jr = j^cos"0, (23) 

to the zonal momentum part of Eq. (^). Here, n is an integer and is set to unity in the 
runs presented in this paper; we have checked that narrower forcing, given by n = 2 or 4, 
leads to qualitatively similar results. Note that the forcing is zonally-symmetric, even though 
the underlying baroclinic process modeled is zonally-asymmetric. As already discussed, the 
exact (and unknown) form of forcing need not be specified and is actually not sought: our 
aim here is an assessment of the evolution of an established superrotating equatorial jet and 
the nature of its interaction with eddies in an idealized setting. Accordingly, the forcing is 
assumed not to respond to the modeled layer itself (i.e., no back- reaction) ; we have in mind 
an "active" modeled layer that overlies a massive "abyssal" interior, which merely acts as 
a source of low angular momentum and friction (e.g., due to convective turbulence). One 
can think of this forcing as arising from some general zonally-asymmetric baroclinic process 
that produce eddy-momentum convergence at low latitude. A small amount of linear drag, 
of the form —au, is added to balance the forcing. 

The resulting flow from Model Zl in Table [2] is depicted in Figure [T71 The resolution of 
the simulation is T106. The iso-g contours in cylindrical projection, centered on the substel- 
lar point. Three time frames are shown: r = {35,40, 100} (Figure [T7k-c. respectively). The 
corresponding zonal jet profiles in each time frames are shown as well. In this simulation, 
the force-dissipation b alance is set so that the maxi mum wind is ~1000 m s~^, similar to the 
simulation reported in lShowman and GuillotI (120021) at the 5.6-8 .1 ba r pressure level. A qui- 
escent initial state is adopted, as in IShowman and GuillotI (120021 ) and ICooper and Showman 
( 120051 ). In the fi gure, the first thing to note i s that the flow is meridionally (north-south) 
symmetric, as in IShowman and GuillotI (120021 ). This is expected given the setup and indi- 
cates a complete dominance of the applied forcing. Particularly clearly visible in Figure [T7k 
are the large Rossby waves that have been excited by the forcing. One can also clearly 
see the superrotation ([m(0 = 0)] > 0) that is produced by the forcing in the correspond- 
ing jet profile. However, as the Rossby waves propagate meridionally, they interact with 
the background flow to drag the mid-latitude mean flow back (in the westward direction), 
simultaneously causing the equatorial mean flow (the source location) to speed up. 

The behavior follows a well known acceleration/deceleration mechanism in atmospheric 
dynamics. The dissipation occurs in the region where the phase speed of the waves match the 
mean flow speed. Since Rossby waves propagate westward, the flow must be eastward for the 
mechanism to operate. Observe in Figure [T7b that, where the wave has dissipated, large-scale 
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eddies have formed. The flow is at this point very similar to that of IShowman and Guillot 



( 120021 ) (see their Figure 6). However, ultimately, the flow is not stable — at least in the 



angular momentum forced equivalent-bar otropic system simulated he re (Figure [T7b). The 



resulting flow is very complex (compared to lShowman and Guillotll2002l ) . but the hemispheric 



symmetry still persists; the increased complexity is due to the higher resolution. One can 
see the strong influence of the propagating waves in Figure [T7b . after which the jet proflle 
does not change very much. Except for the equatorial jet, the mean flow is wholly westward 
from the time shortly before r=100. Clearly, the strength of superrotation depends on the 
pre-existing background flow, given the forcing amplitude. 

Using the equivalent-barotropic model, we have also explored the behavior of super- 
rotating flow states under various initial conditions. Our equihbrated-state flndings are 
summarized in Figure [18], which present results from Models Z2, Zl, and Z3 in Table O 
The flgure illustrates the importance of allowing eddies to be present in the flow fleld. The 
projections for g- fleld maps on the left column are as in Figure [T71 

Figure [TSh illustrates thermally forced (?7 = 0.4) flow without applied angular momentum 
forcing (JF = 0). The jet proflle consists of a mid-latitude westward jet in each hemisphere, 
consistent with spreading of q lines at mid-latitudes by the applied day-night forcing. The 
flow is not very "dynamic" in this case and does not change much over time. In contrast, 
with both angular momentum and day-night forcings applied, the flow, shown in Figure [T8b. 
becomes more complex and dynamic, with the background mean flow interacting strongly 
with radiating planetary (Rossby) waves. Finally, we have in Figure [TSb all three forcings 
(day-night, angular momentum, and initial random stirring) included in the simulation. 
Here, the flow is similar to the turbulent flow cases we have presented throughout this 
paper, with f7 7^ 0. In this case, there are three broad jets. Because the U associated with 
the stirring is 400 m s~^, which is larger than the maximum jet amplitude without the 
stirring (see Figure [TSb). the equatorial jet here is weakly retrograde. Hence, it appears that 
the strength of superrotation is sensitive to the vigor of eddies present. More broadly, the 
absence of vigorous eddies may thus be a general limitation in dynamical simulations. 



6. Conclusion 



In this work, we have extensively explored global, turbulent, adiabatic atmospheric 
dynamics of extrasolar giant planets in circular orbits. Extrasolar giant planets with ec- 
centricity < 0.05 comprise a signiflcant fraction of the currently known population. Several 
robust propertie s emerge from our set of simulations, which broadly conflrm the results in 
Cho et al.l (120031 ). First, there is a strong rotational control on the flow dynamics. This leads 
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to atmospheric flows with a strong zonal component on extrasolar giant planets, independent 
of the thermal forcing amplitude. Hence, the velocity field is not one of a simple diverging 
flow away from the substellar point. This is so even for close-in extrasolar giant planets. This 
behavior is reminiscent of Uranus, which possesses strong zonal winds like the other Solar 
System giant planets even though it is heated at the pole. Second, in trying to understand 
the general circulation of extrasolar giant planet atmospheres, it is clearly useful to divide 
the task along the following two lines: i) understand the zonally-symmetric and zonally 
asymmetric circulation; and, ii) understand the tropical-subtropical and midlatitude-polar 
dynamics. Of course, each includes the interaction with the other two components. Third, 
in the adiabatic equivalent-barotropic equations model, it is possible to have a full range 
of global temperature distributions. They include rotating, oscillating, shifted, and fixed 
day-night distributions. 

Not surprisingly, the spatiotemporal behavior exhibited by the temperature distribution 
depends on the amplitude of the thermal forcing. When thermal forcing is strong, the day- 
night temperature difference is expected to be fixed. In effect, the temperature field decouples 
from the dynamics and can essentially be thought of as "thermal orography" . However, the 
degree of decoupling depends on the strength of the background flow so that, if pre-existing 
zonal winds and eddies are strong, the temperature field is again enslaved to the flow field. 
Therefore, it appears from the adiabatic calculations performed in this work that it would be 
difficult to ascertain the characteristic flow speed independently of the radiative properties 
(e.g., albedo) — and vice versa. Note that because not as much irradiation is absorbed at 
higher altitudes (in the absence of special absorbers), the smaller amplitude forcing situations 
can be considered very roughly applicable to those regions. This would lead to a vertical 
shear in the temperature field and a baroclinic adjustment might be indicated. It would 
be interesting to study such adjustment processes when vertical communication between 
different layers is taken into account. 

In general, the response of the atmospheric motion to both mechanical and radiative 
types of forcing is very complex and difficult to analyze from first principles. Modeling 
such complexity requires a range of approaches — from simple analytical calculations to full 
general circulation modeling. Different approaches are more successful with different pieces of 



the fu ll problem, and th e y are all needed to make good progre s s (e.g . , IShowman and Guillot 



( 2002 ): ICho et all ( 120031 ): iMenou et all (l2003f ): lBurkert et al. I (l2005f ): ICooper and Showman 



(120051 )). In this work, the equivalent-barotropic equations in isentropic coordinate have been 
used to focus on lateral dynamics near the top of the planet's convection zone. The isentropic 
layers above this region do not lie at constant geometric heights. The layers slant slightly 
downward from the substellar point to the poleward direction on extrasolar giant planets; on 
close-in extrasolar giant planets the layers slant downward toward the antistellar direction 
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as well. On the other hand, the isentropic layers below the region (the top of the convection 
zone) slant upward by a very significant amount from the substellar point. In general, air 
parcels move both along and across these surfaces, transporting heat in both horizontal 
and vertical directions. We have taken advantage of the clear association of horizontal and 
vertical motions with adiabatic and diabatic heating, respectively, to perform equivalent- 
barotropic equations simulations in what amounts to "l|-layers" (i.e., a layer with variable 
thickness) . 

Despite the vertical integration, the equivalent-barotropic equations support many of 
the phenomena and types of fluid motion supported by the primitive equations. This in- 
cludes Rossby waves, gravity waves, balanced motions (e.g. geostrophic) , adjustments, and 
barotropic instability. Some of the consequences of baroclinic processes, such as stirring by 
eddies or convection can also be represented and spatially reasonably resolved. The high 
resolution also allows mixing of q and fine-scale tracers and turbulent cascades to small 
scales — all critical features in real atmospheres. Consequences of thermal forcing on eddies 
and mean-fiow can be included through the deflection of lower boundary height. In addi- 
tion, the equivalent-barotropic equations have a consistent set of conservation laws for mass, 
energy, angular momentum, g, potential enstrophy (TYg^), and more exotic quantities like 
pseudomomentum (wave activity). In view of all these properties, we have carefully explored 
the dynamical behavior of the adiabatic equivalent-barotropic equations model in this work. 

In subsequent studies, we will consider more general physical situations important for 
understanding atmospheric dynamics on cxtrasolar giant planets. We will explore the conse- 
quences of relaxing the important adiabatic and barotropic assumptions made in the present 
work, by using models that explicitly include radiative forcing and vertical coupling between 
multiple atmospheric layers. With the recent breakthrough detections of infrared emission 
from several close-in extrasolar giant planets, observational programs should be able to pro- 
vide better constraints in the future on key atmospheric parameters, allowing more accurate 
assessments of the circulation and thermal structures on these planets. 
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Table 1: Global Planetary Parameters 



Planet'^' 




p(3) 
orb 




e(5) 




4"' 


g(8) 




rip 








(days) 


(AU) 




(Mj) 


(m) 


(m 


(rad B -1) 


(m) 


(m 


Jupiter 


1.0 


4,332.6 


5.2 


0.0489 


1.0 


7.1 X lO'' 


23 


1.8 X 10""* 


2 X 10* 


70 


Saturn 


1.0 


10,759.2 


9.58 


0.0565 


0.3 


6 X lO'' 


9 


1.6 X 10""* 


4 X 10* 


400 


Uranus 


1.0 


30,685.4 


19.2 


0.0457 


0.046 


2.6 X lO'' 


9 


(-)l X 10-* 


3.5 X 10* 


300 


Neptune 


1.0 


60,189.0 


30.05 


0.0113 


0.054 


2.5 X lO'' 


11 


9.75 X lO-'* 


3 X 10* 


300 


HD 209458 b 


1.05 


3.5247 


0.045 


0.0 


0.69 


10^ 


8 


2.1 X 10"^ 


7 X 10^ f*"' 





Note. — (1) Giant planets (2) Parent star mass (3) Orbital period (4) Semi-major axis (5) Eccentricity 
(6) Planet mass (7) Planet radius (8) Surface gravity (9) Rotation rate (10) Global pressure scale height (11) 
Global root mean square velocity scale (a) Assuming spin-orbit synchronization (b) Prom global radiative 
equilibrium. 



Table 2: Summary of Simulations Discussed 



Model [/(I) 


To (2) ^(3) 




bands 


Remarks 


(m s~^) 


(K) 




(Rp) 





SI 


70 


130 





.03 


.23 


14 


Jupiter 


S2 


400 


95 





.03 


.64 


5 


Saturn (o-map not shown) 


S3 


300 


60 





.11 


1.1 


3 


Uranus (g'-map not shown) 


S4 


300 


60 





.11 


1.1 


3 


Neptune 


Jl 


70 


130 





.02 


.23 


14 


SI w/ different size initial stirring 


J2 


70 


130* 





.02 


.23 


14 


Jl w/ shallow- water model 


Ul 


70 


260 





.05 


.23 


14 


"Warm" Jupiter 


U2 


400 


1450 





1.2 


1.2 


3 


"Unsynchronized"' HD 209458 b 


U3 


400 


1800 





1.3 


1.2 


3 


"Warm" U2 


U4 


400 


800 





0.9 


1.2 


3 


"Cool" U2 


HI 


100 


1450 


.0 


1.2 


0.6 


5 


HD 209458 b, no thermal contrast 


H2 


100 


1450 


.02 


1.2 


0.6 


5 


HD 209458 b, low thermal contrast 


H3 


100 


1450 


.04 


1.2 


0.6 


5 


HD 209458 b, medium contrast 


H4 


100 


1450 


.20 


1.2 


0.6 


5 


HD 209458 b, high thermal contrast 


H5 


200 


1450 


.04 


1.2 


0.8 


4 


HD 209458 b, low speed jets 


H6 


400 


1450 


.04 


1.2 


1.2 


3 


HD 209458 b, medium speed jets 


H7 


1000 


1450 


.04 


1.2 


1.8 


2 


HD 209458 b, high speed jets 


Zl 





1450 


.04 


1.2 


1.8i 


2 


H6 w/ zero initial vel. &; superrot. 


Z2 





1450 


.04 


1.2 


1.8t 


2 


Z2 w/o superrot. 


Z3 





1450 


.04 


1.2 


1.8t 


2 


Zl w/ eddies 



Note. — (1) Initial global root mean square velocity (2) Global mean temperature at cloudtop for Solar 

System giant planets and fiducial cloudtop for extrasolar giant planets (3) Forcing amplitude (4) Nondimen- 
sional Rossby deformation radius at the pole, based on Hp (5) Nondimensional Rhines length at the equator, 
based on U (6) Number of zonal jets estimated, based on (f ) May eventually erode to 2, if Ln. k, 0.3 (*) 
Based on layer thickness (|) Based on U at quasi-equilibrium. 
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Fig. 1. — Equivalent-barotropic turbulence simulation of Jupiter. Contour maps of the flow 
tracer, potential vorticity g, from six time frames are shown in orthographic view centered 
at the equator. Positive (negative) values are in full (dashed) contours; 40 contour levels 
are shown. Time, in unit of planetary rotation periods (r = 27r/r2), is indicated in the 
upper left corner of each frame. Using only five physical parameters, which are known 
from observations, the initially random turbulent flow self-organizes into one dominated by 
zonal (east-west) bands — alternating, high/low q gradients in the meridional (north-south) 
direction — similar to the actual Jupiter. 

Fig. 2. — Simulation of Figure [1] in stereographic view. The boundary of the disk is the 
equator, and the center of the disk is the north pole. In this projection, the polar region 
occupies a smaller area of the disk compared with the equatorial region, as in Figure [H A 
quasi-steady banded state is reached and is robust. The general flow picture does not change 
even at time, r = 2000, the duration of this run. The interaction of the flow structures 
(vortices and jets) is weak, due to the extremely small Rossby deformation radius, Lti (see 
text for definition). Jets and vortices do not coalesce. 



Fig. 3. — Quasi-steady zonal jet (positive eastward) profile from the simulation of Figures [T] 
and [2] at r = 300. The jets are associated with the bands in Figures [1] and [2l Several 
quahtative features (number, width, and strengths of jets) of the observed profile on Jupiter 
are reproduce d, giving confiden c e in th e model approach. As in the shallow-water model 
calculations of ICho and Polvanil (Il996bl ). however, the sign of the equatorial jet is opposite 
to that observed. This is a common feature of one-layer turbulent models of Solar System 
giant planets without forcing. The observed prograde (eastward) equatorial jet is due to 
mechanism(s) not included in these unforced models. 



Fig. 4. — Detailed comparison with the shallow-water model. Shallow-water simulation (a) 
and corresponding equivalent-barotropic simulation (b) are essentially identical. However, 
the equivalent-barotropic model is superior to the shallow-water model, since crucial param- 
eters (e.g., the equivalent depth He) take on a more physical meaning. Note the presence of 
a large, stable anticyclone (a vortex, which is counter-rotating with respect to fi, which is 
ant i- clockwise in the figure) at mid latitude [at ~315° in (a) and at ~230° longitude in (b) 
at ~40° latitude in both]. This structure is similar to Jupiter's Great Red Spot and is of 
approximately the same size and shape. The structure, along with the bands and jets, also 
appears in Figure [21 which has a different initial condition. Hence, these are robust features 
in our simulations. 



Fig. 5. — Simulation of Neptune. In contrast to Jupiter (e.g., Figures [T] and [2]), there is a large 
band of homogenized potential vorticity around the equator and a strong circumpolar vortex 
centered at each pole (r = 146 frame). A very small number of contours at low latitudes is 
present, compared with the high latitudes. The overall flow pattern is a direct consequence 
of the larger U and smaller Q values (i.e., larger L-ji/Rp and L/^/Rp, respectively), compared 
to those of Jupiter. 
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Fig. 6. — Zonal jet profiles of the Solar System giant planets: (a) Jupiter, (b) Saturn, 
(c) Uranus, and (d) Neptune. The qualitative features of the jets, which are the principle 
structures of the large-scale circulation on the Solar System giant planets, are captured for 
all four planets. Only the observed parameters of Table [1] are used, and no explicit thermal 
forcing is imposed. The width of the jets roughly correspond to L^. On all four planets, 
Ln/Rp ^1/3, which allows the formed jets to be stable over long time. The sign of the 
equatorial jets on Jupiter and Saturn are reversed, compared with observations; but, the 
sign for Uranus and Neptune is as observed. The scale for Jupiter is different from that for 
the others. 
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Fig. 7. — A "warm" Jupiter, with an equilibrium temperature, Tg = 300 K. Contour maps 
of the potential vorticity (g) field at r = 199 are shown in equatorial orthographic (a) and 
polar stereographic (b) views. The physical parameters used are identical to the simulation 
in Figures [1] and [2|, except that is 2.3 times larger (see text for explanation). This case can 
serve as an example of a "Jupiter" which has migrated in to ~1 AU from the central star. 
Additionally, it could apply to a "cool" , unsynchronized extrasolar giant planet with Jupiter 
f2 or to a level 2.3 pressure scale heights deeper than in previous models (compare with 
Figures [1] and [2]) . Due to the increased L-ji value, the bands and jets are less pronounced, 
and the equatorial jet amplitude is reduced. There are also fewer vortices in the polar region. 
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Fig. 8. — Simulation of a HD 209458 b-like, unsynchronized close-in extrasolar giant planet 
in polar stereographic view: early evolution. The rotation rate for this planet is that of 
HD 209458 b (i.e., 8.4 times slower than that of Jupiter). The global average root mean 
square velocity f7 for HD 209458 b is not known and can range plausibly between ~10^ m s~^ 
to ~10^ m s~^. In this simulation, U = 400 m s~^, roughly the maximum value of all the 
planets in the Solar System. We have varied this parameter in our study and summarize 
the findings in §5. The salient feature here is that the vortices and jets in the flow are much 
more dynamic compared with those in Figure Ej and even with those in Figure [5l This is 
due to the larger Ln value, which is ~-Rp in this case. 
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Fig. 9. — Same simulation and view of Figure [HI evolution post equilibrium. By r = 70 a 
broad band of smoothly varying q forms outside a robust polar vortex, which has formed 
through continuous mergers. The flow is more similar to that on Neptune than on Jupiter 
(compare with Figures [2] and [5]) . As in Neptune, a polar vortex forms in each hemisphere, 
but the vortex is not centered on the pole. Robust, off-pole vortices are general features 
in simulations of the circulation of HD 209458 b. Another difference, consistent with less 
homogenized equatorial band, is the planetary (Rossby) wave breaking that appears more 
pronounced at mid-latitudes. 
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Fig. 10. — Quasi-steady jet profiles corresponding to the simulation of Figures [S] and M (a) 
and the same simulation with a different U, which is 100 m s'^ (b). The low number of jets 
is also a robust feature of our synchronized HD 209458 b simulations. The precise profile of 
the jets depends on the value of U, and the profile is slightly time-varying in a given run. 
With a value of U characteristic of most Solar System giant planets, three or four broad 
jets form. When U ^ 100 m s~^ more jets may form. Occasionally a two-jet profile forms 
for U ^ 400 m s~^ and also for U = (see Figure [T7|) . Note that in (b) the equatorial jet 
is prograde, demonstrating that a retrograde equatorial jet is not a necessary outcome of 
single-layer simulations. That outcome requires a small L-ji value. 
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Fig. 11. — Unsynchronized HD 209458 b simulations with different equilibrium temperatures, 
Tq: (a) 1800 K and (b) 800 K. The two temperatures are chosen to roughly bracket a variety 
of physical situations (e.g., distance, stellar flux, and albedo). For fixed flux and albedo, the 
higher Tq value can also represent an HD 209458 b-like planet which is ~.005 AU closer to 
the host star than HD 209458 b is to its star. Similarly, the lower Tq value can represent 
an HD 209458 b-like planet that is ~0.15 AU further away. The higher Tq value leads to a 
more smoothly varying low- to mid-latitude potential vorticity distribution, while the lower 
Tq value leads to much more planetary (Rossby) wave activity, along with a more centered 
polar vortex. The basic three-jet profile is not changed in either case. 
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Fig. 12. — HD 209458 b atmospheric flow for different forcing amplitudes, rj: (a) 0., (b) 
.02, (c) .04, and (d) .20. All other parameters in (a)-(d) are identical; in particular, U = 
100 m s"^. Formally, rj can range from to 1 — i.e., from no asymmetric forcing to full forcing, 
respectively. The latter corresponds to a temperature deviation, 6T — 0{To), which is not 
physically realistic on large scales. The primary effect of increasing the forcing asymmetry 
is to sharpen the midlatitudc jet (a-c). When the forcing is large, it induces large zonal 
asymmetry in the jet, destabilizing it and causing waves to break at high latitudes (d). A 
simulation with 77 = 0.4 shows a flow that is basically a more violent version of (d). 



Fig. 13. — Temperature distribution, Tq = To(A,(/)), of the runs in Figure [T2] with rj: (a) 
0., (b) .02, (c) .04, and (d) .20. For a fixed U, Tq is strongly dependent on rj. With no or 
low-amphtude forcing, Tq essentially tracks the fiow structure (a). As the forcing amplitude 
increases, however, there is a competition between the fiow-induced and the forcing- induced 
temperature patterns: temperature anomalies either rotate quasi-steadily about the pole (b) 
or oscillate aperiodically about a point near the substellar point (c). For high rj values, the 
forcing dominates, and Tq is essentially the day-night distribution centered on the substellar 
point (d). In each case, a cyclone/anticyclone pair at high latitudes is present (i.e., inde- 
pendent of the amplitude) and is a robust feature of our simulations (see Figures [12] and 
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Fig. 14. — HD 209458 b circulation for different global root mean square velocities, U : (a) 
100 m s~^, (b) 200 m s"^, (c) 400 m s~^, and (d) 1000 m s~^. The primary effect of larger U 
values is to strengthen and enlarge the polar vortex at high latitudes. At low to mid latitudes, 
breaking of planetary (Rossby) waves appear to be reduced when U is larger. However, in 
reality, the breaking also increases, as shown in Figure [151 The increase in vortex size is due 
to the initially more energetic vortices, which coalesce to form the final vortex. The vortex 
in (d) is large enough to nearly cap the planet, poleward of 60 degrees latitude, if it were 
centered at the pole. 
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Fig. 15. — Flow outside the polar vortex and global jet profile (for two different U values). 
Three same time frames of the simulation in Figure UMi are shown in (a-c). There are 50 
contours in each map. Note the dominance of the cyclonic polar vortex (a) and the presence 
of a weaker, anticyclonic vortex, which is revealed when a smaller contour range is used 

(b) . The anticyclone forms from high-amplitude breaking flow outside the cyclonic vortex 

(c) . The cyclone is associated with a coherent cold region, and the anticyclone is associated 
with a coherent warm region — as shown in Figure [ini Zonal jet profiles in two runs with 
different U values: (d) 100 m s~^ and (e) 1000 m s^^. The number and width of the jets is 
qualitatively consistent with the Rhines scale, L^. 



Fig. 16. — HD 209458 b atmospheric temperature distribution Tq for different global root 
mean square velocity U : (a) 100 m s~^, (b) 200 m s^^, (c) 400 m s~^, and (d) 1000 m s~^. The 
same simulations as in Figure [T3] are shown. With forcing amplitude rj fixed, Tq is strongly 
dependent on If. At low If values, the forcing dominates and Tq is essentially the day-night 
distribution centered on the substellar point (a). As U increases, competition between the 
forcing-induced and the flow-induced temperature patterns develops: temperature anomalies 
oscillate aperiodically about a point near the substellar point (b) or rotate quasi-steadily 
about the pole (c). At large If values, Tq essentially tracks the flow structure on a several- 
days timescale (d). 
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Fig. 17. — Simulation starting at rest. The model is augmented to include the baroclinic 
effect of large-amplitude, zonally-asymmetric, equatorial forcing. Potential vorticity (q) con- 
tours in cylindrical projection, centered at the substellar point, are shown for three time 
frames (a-c). Corresponding zonal jet profiles are shown to the right of each q-map. The 
forcing excites planetary (Rossby) waves (a). Angular momentum is transported as propa- 
gating planetary waves interact with the background flow (b). Eventually, the flow becomes 
very complex, at the high resolution of this simulation (c). However, the north-south sym- 
metry is maintained throughout, indicating the dominance of the forcing. 
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Fig. 18. — Summary of results with quiescent initial condition. The plotting is identical 
to Figure [171 Late-time results in three simulations with different forcing conditions are 
shown: (a) day-night thermal forcing only, (b) with angular momentum forcing added, and 
(c) with initial random stirring added (i.e. atmosphere no longer at rest initially). With 
only day-night forcing applied, the flow is laminar and consistent with q lines distorted by 
the thermal forcing. Finally, with angular momentum forcing also added, the flow becomes 
more complex and the atmosphere superrotates. With the addition of initial random stirring, 
the flow becomes "eddy driven" , if the stirring is energetic. The equatorial jet is weakly 
retrograde in this case. 



